# Project Information---
#  Project: Educating for Inclusion: Diversity education programs can reduce prejudice towards outgroups in Israel
#  Journal: PNAS
#  Last updated: Tue Mar  7 15:42:43 2023
#  Purpose: Estimate all results for paper -- Study II
#  Outputs: All figures and tables associated with Study 2 in paper and appendix
#  Machine: Chagai's Macbook Pro



#Load Packages------
library("tidyverse")
library("pmdplyr")
library("DataCombine")
library("dplyr")
library("readxl")
library("lubridate")
library("estimatr")
library("ltm")
library("effectsize")
library("ggthemes")
library("stargazer")
library("fwildclusterboot")
library("RItools")
library("xtable")
library("multiwayvcov")
library("lmtest")
library("RColorBrewer")
library("interflex")
library("ri2")
library("exr")
library("forcats")

set.seed(999)

#Read clean data------
setwd(paste(getwd()))
full_data <- 
 readRDS("study2_data.rds")


# Results in Main Text-------
##Figure 4: Main ATEs------
###Thermometor------
therm_index <- 
 full_data %>% 
 lm_lin(standardize(therm_index_post) ~ 
         treatment,
        covariates = ~  therm_index_pre+
         block_id +male + 
         na_therm_arab_pre + na_therm_blind_pre + 
         na_therm_imgrnt_pre + na_therm_uo_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
therm_index_out <-
 tidy(therm_index) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Target = "Index",
        Outcome = "Thermometers",
        Wave = "1-2 Weeks")


therm_index2 <- 
 full_data %>% 
 lm_lin(standardize(therm_index_post2) ~ 
         treatment,
        covariates = ~  therm_index_pre+
         block_id +male + 
         na_therm_arab_pre + na_therm_blind_pre + 
         na_therm_imgrnt_pre + na_therm_uo_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
therm_index_out2 <-
 tidy(therm_index2) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Target = "Index",
        Outcome = "Thermometers",
        Wave = "8-13 Weeks")

###Contact------
cntct_index <- 
 full_data %>% 
 lm_lin(standardize(cntct_indx_post) ~ 
         treatment,
        covariates = ~  cntct_indx_pre+
         block_id +male +
         na_cntct_1_blind_pre + na_cntct_2_blind_pre +
         na_cntct_3_blind_pre + na_cntct_4_blind_pre +
         na_cntct_1_arab_pre + na_cntct_2_arab_pre + 
         na_cntct_3_arab_pre + na_cntct_4_arab_pre +
         na_cntct_1_imgrnt_pre + na_cntct_2_imgrnt_pre +
         na_cntct_3_imgrnt_pre + na_cntct_4_imgrnt_pre +
         na_cntct_1_uo_pre + na_cntct_2_uo_pre +
         na_cntct_3_uo_pre + na_cntct_4_uo_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
cntct_index_out <-
 tidy(cntct_index) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Outcome = "Contact \nIntentions",
        Target = "Index",
        Wave = "1-2 Weeks")


cntct_index2 <- 
 full_data %>% 
 lm_lin(standardize(cntct_indx_post2) ~ 
         treatment,
        covariates = ~  cntct_indx_pre+
         block_id +male+ 
         na_cntct_1_blind_pre + na_cntct_2_blind_pre +
         na_cntct_3_blind_pre + na_cntct_4_blind_pre +
         na_cntct_1_arab_pre + na_cntct_2_arab_pre + 
         na_cntct_3_arab_pre + na_cntct_4_arab_pre +
         na_cntct_1_imgrnt_pre + na_cntct_2_imgrnt_pre +
         na_cntct_3_imgrnt_pre + na_cntct_4_imgrnt_pre +
         na_cntct_1_uo_pre + na_cntct_2_uo_pre +
         na_cntct_3_uo_pre + na_cntct_4_uo_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
cntct_index_out2 <-
 tidy(cntct_index2) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Outcome = "Contact \nIntentions",
        Target = "Index",
        Wave = "8-13 Weeks")

###Contact Behavior------
cntct_bhvr <- 
 full_data %>% 
 lm_lin(standardize(contact_project) ~ 
         treatment,
        covariates = ~  cntct_indx_pre+ therm_index_pre + div_scale_pre+
         block_id +male +
         na_cntct_1_blind_pre + na_cntct_2_blind_pre +
         na_cntct_3_blind_pre + na_cntct_4_blind_pre +
         na_cntct_1_arab_pre + na_cntct_2_arab_pre + 
         na_cntct_3_arab_pre + na_cntct_4_arab_pre +
         na_cntct_1_imgrnt_pre + na_cntct_2_imgrnt_pre +
         na_cntct_3_imgrnt_pre + na_cntct_4_imgrnt_pre +
         na_cntct_1_uo_pre + na_cntct_2_uo_pre +
         na_cntct_3_uo_pre + na_cntct_4_uo_pre+
         na_therm_arab_pre + na_therm_blind_pre + 
         na_therm_imgrnt_pre + na_therm_uo_pre+
         na_div1_pre + na_div2_pre + na_div3_pre +
         na_div4_pre + na_div5_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
cntct_bhvr_out <-
 tidy(cntct_bhvr) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Outcome = "Register \nContact",
        Target = "Behavioral \nMeasure",
        Wave = "1-2 Weeks")



cntct_bhvr2 <- 
 full_data %>% 
 lm_lin(standardize(contact_project_post2) ~ 
         treatment,
        covariates = ~  cntct_indx_pre+ therm_index_pre + div_scale_pre+ 
         block_id +male + 
         na_cntct_1_blind_pre + na_cntct_2_blind_pre +
         na_cntct_3_blind_pre + na_cntct_4_blind_pre +
         na_cntct_1_arab_pre + na_cntct_2_arab_pre + 
         na_cntct_3_arab_pre + na_cntct_4_arab_pre +
         na_cntct_1_imgrnt_pre + na_cntct_2_imgrnt_pre +
         na_cntct_3_imgrnt_pre + na_cntct_4_imgrnt_pre +
         na_cntct_1_uo_pre + na_cntct_2_uo_pre +
         na_cntct_3_uo_pre + na_cntct_4_uo_pre+
         na_therm_arab_pre + na_therm_blind_pre + 
         na_therm_imgrnt_pre + na_therm_uo_pre+
         na_div1_pre + na_div2_pre + na_div3_pre +
         na_div4_pre + na_div5_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
cntct_bhvr_out2 <-
 tidy(cntct_bhvr2) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Outcome = "Register \nContact",
        Target = "Behavioral \nMeasure",
        Wave = "8-13 Weeks")




###Diversity------
div_index <- 
 full_data %>% 
 lm_lin(standardize(div_scale_post) ~ 
         treatment,
        covariates = ~  div_scale_pre+
         block_id +male + 
         na_div1_pre + na_div2_pre + na_div3_pre +
         na_div4_pre + na_div5_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
div_index_out <-
 tidy(div_index) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Outcome = "Diversity \nAttitudes",
        Target = "Index",
        Wave = "1-2 Weeks")

div_index2 <- 
 full_data %>% 
 lm_lin(standardize(div_scale_post2) ~ 
         treatment,
        covariates = ~  div_scale_pre+
         block_id +male + 
         na_div1_pre + na_div2_pre + na_div3_pre +
         na_div4_pre + na_div5_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
div_index_out2 <-
 tidy(div_index2) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Outcome = "Diversity \nAttitudes",
        Target = "Index",
        Wave = "8-13 Weeks")



### Diversity Behavior------

div_bhvr2 <- 
 full_data %>% 
 filter(.,
        school != "Democrati") %>% 
 lm_lin(standardize(wrist_band) ~ 
         treatment,
        covariates = ~  cntct_indx_pre+ therm_index_pre + div_scale_pre+ 
         block_id + male +
         na_cntct_1_blind_pre + na_cntct_2_blind_pre +
         na_cntct_3_blind_pre + na_cntct_4_blind_pre +
         na_cntct_1_arab_pre + na_cntct_2_arab_pre + 
         na_cntct_3_arab_pre + na_cntct_4_arab_pre +
         na_cntct_1_imgrnt_pre + na_cntct_2_imgrnt_pre +
         na_cntct_3_imgrnt_pre + na_cntct_4_imgrnt_pre +
         na_cntct_1_uo_pre + na_cntct_2_uo_pre +
         na_cntct_3_uo_pre + na_cntct_4_uo_pre+
         na_therm_arab_pre + na_therm_blind_pre + 
         na_therm_imgrnt_pre + na_therm_uo_pre+
         na_div1_pre + na_div2_pre + na_div3_pre +
         na_div4_pre + na_div5_pre,
       clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
div_bhvr_out2 <-
 tidy(div_bhvr2) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Outcome = "Diversity \nWristband",
        Target = "Behavioral \nMeasure",
        Wave = "8-13 Weeks")


#### Plot main results------

summary_study2 <- rbind(therm_index_out, therm_index_out2,
                        cntct_index_out, cntct_index_out2,
                        cntct_bhvr_out, cntct_bhvr_out2,
                        div_index_out, div_index_out2,
                        div_bhvr_out2) %>% 
 mutate(.,
        se = paste0("(", round(std.error,2),")"),
        pese = paste(round(estimate,2), se))

# set order
summary_study2$Outcome <- 
 factor(summary_study2$Outcome, 
        c("Diversity \nWristband", "Register \nContact", 
          "Diversity \nAttitudes",
          "Contact \nIntentions",
          "Thermometers"))

# Add annotation
annotation <- data.frame(
 x = c(-0.1,-0.1),
 y = c(5.2,4.8),
 label = c("1-2 Weeks", "8-13 Weeks")
)

ggplot(summary_study2, aes(x=estimate, y=Outcome,
                           colour = forcats::fct_rev(Wave))) + 
 geom_point(position = position_dodge(width = 0.6))+
 geom_errorbar(aes(xmin=conf.low, xmax=conf.high),
               position = position_dodge(width = 0.6),
               width=.1) +
 geom_vline(xintercept = 0,
            color = "black",
            linetype = "dashed",
            size = 0.4) +
 xlim(-0.18,.7)+
 #        geom_label()+
 geom_text(aes(label= pese, color = forcats::fct_rev(Wave)), 
           position = position_dodge(width = .6),
           family = "Times",
           hjust = -1, 
           size = 3,
           show.legend = FALSE) +
 geom_label(data=annotation, aes(x=x, y=y, label=label),
            color=c("dodgerblue4", "firebrick3"), 
            size=2 , angle=45, fontface="bold")+
 labs(x = "",
      y = "",
      color = "")+
 scale_color_manual(values = c("firebrick3", "dodgerblue4"))+
 # theme_tufte()+
 guides(color = guide_legend(reverse = TRUE))+
 theme(panel.spacing = unit(1.5, "lines"),
       text = element_text(size = 15, family = "Times"),
       # axis.text.y = element_text(size=14),
       legend.position = "none",
       panel.grid.major = element_blank(),
       panel.background = element_blank(),
       legend.key=element_blank())



## Figure 5: Potential Mechanisms------

###Perspective Taking-----
pt_index <- 
 full_data %>% 
 lm_lin(standardize(pt_index_post) ~ 
         treatment,
        covariates = ~  pt_index_pre+
         block_id +male +
         na_pt_arab_pre + na_pt_blind_pre +
         na_pt_imgrnt_pre + na_pt_uo_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
pt_index_out <-
 tidy(pt_index) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Target = "Index",
        Outcome = "Perspective Taking",
        Wave = "1-2 Weeks")


pt_index2 <- 
 full_data %>% 
 lm_lin(standardize(pt_index_post2) ~ 
         treatment,
        covariates = ~  pt_index_pre+
         block_id +male +
         na_pt_arab_pre + na_pt_blind_pre +
         na_pt_imgrnt_pre + na_pt_uo_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
pt_index_out2 <-
 tidy(pt_index2) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Target = "Index",
        Outcome = "Perspective Taking",
        Wave = "8-13 Weeks")

### Heterogeneity------

het_index <- 
 full_data %>% 
 lm_lin(standardize(het_index_post) ~ 
         treatment,
        covariates = ~  het_index_pre+
         block_id +male +
         na_htr_arab_pre + na_htr_blind_pre +
         na_htr_imgrnt_pre + na_htr_uo_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
het_index_out <-
 tidy(het_index) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Target = "Index",
        Outcome = "Group Heterogeneity",
        Wave = "1-2 Weeks")



het_index2 <- 
 full_data %>% 
 lm_lin(standardize(het_index_post2) ~ 
         treatment,
        covariates = ~  het_index_pre+
         block_id +male +
         na_htr_arab_pre + na_htr_blind_pre +
         na_htr_imgrnt_pre + na_htr_uo_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
het_index_out2 <-
 tidy(het_index2) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Target = "Index",
        Outcome = "Group Heterogeneity",
        Wave = "8-13 Weeks")



### Similarity------

sim_index <- 
 full_data %>% 
 lm_lin(standardize(similarity_index_post) ~ 
         treatment,
        covariates = ~  similarity_index_pre+
         block_id +male +
         na_smlr_blind_pre + na_smlr_imgrnt_pre +
         na_smlr_arb_pre + na_smlr_uo_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
sim_index_out <-
 tidy(sim_index) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Target = "Index",
        Outcome = "Group Similarity",
        Wave = "1-2 Weeks")


sim_index2 <- 
 full_data %>% 
 lm_lin(standardize(similarity_index_post2) ~ 
         treatment,
        covariates = ~  similarity_index_pre+
         block_id +male +
         na_smlr_blind_pre + na_smlr_imgrnt_pre +
         na_smlr_arb_pre + na_smlr_uo_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
sim_index_out2 <-
 tidy(sim_index2) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Target = "Index",
        Outcome = "Group Similarity",
        Wave = "8-13 Weeks")



####Create Mechanism Plot-----

mech <- rbind(pt_index_out, pt_index_out2,
               sim_index_out, sim_index_out2,
               het_index_out, het_index_out2) %>% 
 mutate(.,
        se = paste0("(", round(std.error,2),")"),
        pese = paste(round(estimate,2), se))

mech$Outcome <- factor(mech$Outcome,
                       c("Group Heterogeneity", "Group Similarity", "Perspective Taking")) 

annotation_mech <- data.frame(
 x = c(-0.1,-0.1),
 y = c(3.15,2.85),
 label = c("1-2 Weeks", "8-13 Weeks")
)

ggplot(mech, aes(x=estimate, y=Outcome,
                 colour = forcats::fct_rev(Wave))) + 
 geom_point(position = position_dodge(width = 0.6))+
 geom_errorbar(aes(xmin=conf.low, xmax=conf.high), 
               width=.1,
               position = position_dodge(width = 0.6)) +
 geom_vline(xintercept = 0,
            color = "black",
            linetype = "dashed",
            size = 0.4) +
 xlim(-.18,.55)+
 geom_text(aes(label= pese, color = forcats::fct_rev(Wave)), 
           position = position_dodge(width = .6),
           hjust = -1.2, 
           size = 3,
           family = "Times",
           show.legend = FALSE) +
 geom_label(data=annotation_mech, aes(x=x, y=y, label=label),
            color=c("dodgerblue4", "firebrick3"), 
            size=2 , angle=45, fontface="bold")+
 labs(x = "",
      y = "",
      color = "")+
 scale_color_manual(values = c("firebrick3","dodgerblue4"))+
 guides(color = guide_legend(reverse = TRUE))+
 theme(panel.spacing = unit(1.5, "lines"),
       text = element_text(size = 14, family = "Times"),
       legend.position = "none",
       panel.grid.major = element_blank(),
       panel.background = element_blank(),
       legend.key=element_blank())





#Appendix-------

##Figure A1 panels c/d: Levels of pre-treatment thermometor and similarity ----
###Therm
# select relevant measures
therm_data <- full_data %>% 
  dplyr::select(.,
                therm_blind_pre, therm_arab_pre,
                therm_uo_pre, therm_imgrnt_pre)

# pivot data
therm_piv <-
  therm_data %>% 
  pivot_longer(therm_blind_pre:therm_imgrnt_pre, 
               names_to = "group", 
               values_to = "value") 
# get mean & SD
therm_estimates <- 
  therm_piv %>%
  group_by(group) %>%
  summarize(Mean_therm = mean(value, na.rm=TRUE),
            Sd_therm = sd(value, na.rm=TRUE),
            n_therm = n()) %>% 
  mutate(se_therm = Sd_therm / sqrt(n_therm),
         lower_ci_rust = Mean_therm - qt(1 - (0.05 / 2), n_therm- 1) * se_therm,
         upper_ci_rust = Mean_therm + qt(1 - (0.05 / 2), n_therm - 1) * se_therm) %>% 
  mutate(.,
         group = case_when(
           group == "therm_blind_pre" ~ "Visually \nImpaired",
           group == "therm_arab_pre" ~ "Arabs",
           group == "therm_uo_pre" ~ "Ultra-Orthodox",
           group == "therm_imgrnt_pre" ~ "Immigrants"
         ))


ggplot(therm_estimates, aes(x=Mean_therm, y=group, colour="dodgerblue4")) + 
  geom_errorbar(aes(xmin=lower_ci_rust, xmax=upper_ci_rust), width=.1) +
  geom_point()+
  labs(x = "Thermometer",
       y = "")+
  geom_vline(xintercept = mean(therm_piv$value, na.rm = T),
             color = "black",
             linetype = "dashed",
             size = 0.4)+
  xlim(0,100)+
  scale_color_manual(values = c("dodgerblue4"))+
  theme(text = element_text(size = 20, family = "Times"),
        legend.position = "none",
        panel.grid.major = element_blank(),
        panel.background = element_blank())



####Similarity
# select relevant measures
sim_data <- full_data %>% 
  dplyr::select(.,
                smlr_arb_pre, smlr_uo_pre,
                smlr_blind_pre, smlr_imgrnt_pre)

# pivot data
sim_piv <-
  sim_data %>% 
  pivot_longer(smlr_arb_pre:smlr_imgrnt_pre, 
               names_to = "group", 
               values_to = "value") 
# get mean & SD
sim_estimates <- 
  sim_piv %>%
  group_by(group) %>%
  summarize(Mean_sim = mean(value, na.rm=TRUE),
            Sd_sim = sd(value, na.rm=TRUE),
            n_sim = n()) %>% 
  mutate(se_sim = Sd_sim / sqrt(n_sim),
         lower_ci_rust = Mean_sim - qt(1 - (0.05 / 2), n_sim- 1) * se_sim,
         upper_ci_rust = Mean_sim + qt(1 - (0.05 / 2), n_sim - 1) * se_sim) %>% 
  mutate(.,
         group = case_when(
           group == "smlr_blind_pre" ~ "Visually \nImpaired",
           group == "smlr_arb_pre" ~ "Arab",
           group == "smlr_uo_pre" ~ "Ultra-Orthodox",
           group == "smlr_imgrnt_pre" ~ "Immigrant"
         ))


ggplot(sim_estimates, aes(x=Mean_sim, y=group, colour="dodgerblue4")) + 
  geom_errorbar(aes(xmin=lower_ci_rust, xmax=upper_ci_rust), width=.1) +
  geom_point()+
  labs(x = "Similarity Index",
       y = "")+
  scale_color_manual(values = c("dodgerblue4"))+
  geom_vline(xintercept = mean(sim_piv$value, na.rm = T),
             color = "black",
             linetype = "dashed",
             size = 0.4)+
  xlim(1,5)+
  theme(text = element_text(size = 20, family = "Times"),
        legend.position = "none",
        panel.grid.major = element_blank(),
        panel.background = element_blank())


##Figure A3: Study treatment recall-----
full_data <-
  full_data %>% 
  mutate(.,
         lesson1_attendance = ifelse(is.na(lesson1_attendance),0,lesson1_attendance),
         lesson1_attendance = ifelse(is.na(lesson2_attendance),0,lesson2_attendance),
         lesson1_attendance = ifelse(is.na(lesson3_attendance),0,lesson3_attendance),
         lesson1_attendance = ifelse(is.na(lesson4_attendance),0,lesson4_attendance),
         number_vids = lesson1_attendance+ lesson2_attendance+
           lesson3_attendance + lesson4_attendance)

full_data %>% 
  filter(.,
         treatment==1) %>% 
  ggplot(., aes(x=number_vids)) + 
  geom_histogram(fill = "dodgerblue4")+
  labs(x = "Number of Videos Students Recall Watching in Class",
       y = "Count")+
  facet_wrap(~class_id, ncol = 5)+
  theme(text = element_text(size = 14, family = "Times"),
        legend.position = "none",
        panel.grid.major = element_blank(),
        panel.background = element_blank())


##Figure A5: Plot days since intervention-----
days_since <- rbind(cbind(as.numeric(full_data$days_since_intervention), "First Wave"),
                    cbind(as.numeric(full_data$days_since_intervention_post2), "Second Wave")) %>% 
  as.data.frame() %>% 
  dplyr::rename(.,
                Days = V1,
                Wave = V2) %>% 
  mutate(.,
         Days = as.numeric(Days)) %>% 
  filter(.,
         !is.na(Days))

mu <- plyr::ddply(days_since, "Wave", summarise, 
                  grp.mean=mean(Days, na.rm=T))

ggplot(days_since, aes(x=Days, fill=Wave)) +
  geom_histogram(position="dodge",
                 alpha = .6)+
  geom_vline(data=mu, aes(xintercept=grp.mean, color=Wave),
             linetype="dashed")+
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10))+
  labs(x= "Days Since Treatment",
       y = "Count",
       color = "",
       fill = "")+
  scale_fill_manual(values = c("dodgerblue4", "firebrick4"))+
  scale_color_manual(values = c("dodgerblue3", "firebrick3"))+
  theme(panel.spacing = unit(1.5, "lines"),
        text = element_text(size = 14, family = "Times"),
        # axis.text.y = element_text(size=14),
        legend.position = "bottom",
        panel.grid.major = element_blank(),
        panel.background = element_blank(),
        legend.key=element_blank())


##Figure A7 panels a/b/c: Correlation of Prejudice Over Time-----
full_data %>% 
  filter(.,
         treatment==0) %>% 
  ggplot(., aes(therm_index_pre, therm_index_post))+
  geom_point(color = "dodgerblue4")+
  geom_jitter(color = "dodgerblue4")+
  geom_smooth(method='lm', colour = "firebrick4")+
  labs(x="Pre-Treatment Thermometor Index",
       y="Post-Treatment Themometor Index")+
  theme(panel.spacing = unit(1.5, "lines"),
        text = element_text(size = 12, family = "Times"),
        # axis.text.y = element_text(size=14),
        legend.position = "bottom",
        panel.grid.major = element_blank(),
        panel.background = element_blank(),
        legend.key=element_blank())


full_data %>% 
  filter(.,
         treatment==0) %>% 
  ggplot(., aes(cntct_indx_pre, cntct_indx_post))+
  geom_point(color = "dodgerblue4")+
  geom_jitter(color = "dodgerblue4")+
  geom_smooth(method='lm', colour = "firebrick4")+
  labs(x="Pre-Treatment Contact Index",
       y="Post-Treatment Contact Index")+
  theme(panel.spacing = unit(1.5, "lines"),
        text = element_text(size = 12, family = "Times"),
        # axis.text.y = element_text(size=14),
        legend.position = "bottom",
        panel.grid.major = element_blank(),
        panel.background = element_blank(),
        legend.key=element_blank())


full_data %>% 
  filter(.,
         treatment==0) %>% 
  ggplot(., aes(div_scale_pre, div_scale_post))+
  geom_point(color = "dodgerblue4")+
  geom_jitter(color = "dodgerblue4")+
  geom_smooth(method='lm', colour = "firebrick4")+
  labs(x="Pre-Treatment Diversity Index",
       y="Post-Treatment Diversity Index")+
  theme(panel.spacing = unit(1.5, "lines"),
        text = element_text(size = 12, family = "Times"),
        # axis.text.y = element_text(size=14),
        legend.position = "bottom",
        panel.grid.major = element_blank(),
        panel.background = element_blank(),
        legend.key=element_blank())

##Figure A8: Correlation of attitudinal and behavioral prejudice-----
cntct_int <- 
  full_data %>% 
  lm_robust(contact_project~standardize(cntct_indx_post),
            se_type = "stata",
            data = .) %>% 
  tidy(., conf.int = T) %>% 
  filter(.,
         term == "standardize(cntct_indx_post)") %>% 
  mutate(.,
         term = "Contact \nIndex",
         outcome = "Contact Registration")

cntct_therm<-
  full_data %>% 
  lm_robust(contact_project~standardize(therm_index_post),
            se_type = "stata",
            data = .) %>% 
  tidy(., conf.int = T) %>% 
  filter(.,
         term == "standardize(therm_index_post)") %>% 
  mutate(.,
         term = "Thermometer \nIndex",
         outcome = "Contact Registration")

cntct_div <- 
  full_data %>% 
  lm_robust(contact_project~standardize(div_scale_post),
            se_type = "stata",
            data = .)  %>% 
  tidy(., conf.int = T) %>% 
  filter(.,
         term == "standardize(div_scale_post)") %>% 
  mutate(.,
         term = "Diversity \nIndex",
         outcome = "Contact Registration")


### Wrist Band------
wb_int <-
  full_data %>% 
  lm_robust(wrist_band~standardize(cntct_indx_post2),
            se_type = "stata",
            data = .) %>% 
  tidy(., conf.int = T) %>% 
  filter(.,
         term == "standardize(cntct_indx_post2)") %>% 
  mutate(.,
         term = "Contact \nIndex",
         outcome = "Wristband Takeup")

wb_div <-
  full_data %>% 
  lm_robust(wrist_band~standardize(div_scale_post2),
            se_type = "stata",
            data = .) %>% 
  tidy(., conf.int = T) %>% 
  filter(.,
         term == "standardize(div_scale_post2)") %>% 
  mutate(.,
         term = "Diversity \nIndex",
         outcome = "Wristband Takeup")

wb_therm <-
  full_data %>% 
  lm_robust(wrist_band~standardize(therm_index_post2),
            se_type = "stata",
            data = .) %>% 
  tidy(., conf.int=T) %>% 
  filter(.,
         term == "standardize(therm_index_post2)") %>% 
  mutate(.,
         term = "Thermometer \nIndex",
         outcome = "Wristband Takeup")


###Plot Behavior Attitude Correlation----

corr_models <- 
  rbind(cntct_int, cntct_div, cntct_therm,
        wb_int, wb_div, wb_therm)


ggplot(corr_models, aes(x=estimate, y=term, color = "dodgerblue4")) + 
  geom_point(position = position_dodge(width = 0.6), show.legend = FALSE)+
  geom_errorbar(aes(xmin=conf.low, xmax=conf.high), 
                width=.1,
                position = position_dodge(width = 0.6),
                show.legend = FALSE) +
  geom_vline(xintercept = 0,
             color = "black",
             linetype = "dashed",
             size = 0.4) +
  xlim(-.3,.3)+
  # scale_y_discrete(limits = att_order)+
  labs(x = "",
       y = "",
       color = "")+
  scale_color_manual(values = c("dodgerblue4"))+
  facet_grid(cols = vars(outcome))+
  # theme_tufte()+
  #guides(color = guide_legend(reverse = TRUE))+
  theme(panel.spacing = unit(1.5, "lines"),
        text = element_text(size = 12, family = "Times"),
        # axis.text.y = element_text(size=14),
        legend.position = "bottom",
        panel.grid.major = element_blank(),
        panel.background = element_blank(),
        legend.key=element_blank())



##Figure A9: Distribution of behavioral measures------

full_data %>% 
  filter(.,
         treatment!=1) %>% 
  ggplot(., aes(x=contact_project)) + 
  geom_histogram(aes(alpha =0.3),fill = "dodgerblue4")+
  labs(x = "Contact Registration (Control Group First Post Treatment)",
       y = "Count")+
  geom_vline(xintercept = mean(full_data$contact_project[full_data$treatment!=1], na.rm = T),
             color = "dodgerblue4",
             linetype = "dashed",
             size = 0.4)+
  theme(text = element_text(size = 14, family = "Times"),
        legend.position = "none",
        panel.grid.major = element_blank(),
        panel.background = element_blank())




full_data %>% 
  filter(.,
         treatment!=1) %>% 
  ggplot(., aes(x=wrist_band)) + 
  geom_histogram(aes(alpha =0.3),fill = "dodgerblue4")+
  labs(x = "Wristband Take-up (Control Group Second Post Treatment)",
       y = "Count")+
  geom_vline(xintercept = mean(full_data$wrist_band[full_data$treatment!=1], na.rm = T),
             color = "dodgerblue4",
             linetype = "dashed",
             size = 0.4)+
  theme(text = element_text(size = 14, family = "Times"),
        legend.position = "none",
        panel.grid.major = element_blank(),
        panel.background = element_blank())


##Figure 10: Open ended responses about research-related activities in school----
# read in relevant data
open_end <- 
  read.csv("study2_open_end_coded.csv")
open_end <-
  open_end %>% 
  mutate(.,
         coding_with_na = ifelse(coding == "", "No Response", coding),
         coding = ifelse(coding == "",NA, coding),
         response = ifelse(is.na(coding),"No", "Yes"))


open_end %>% 
  filter(.,
         !is.na(coding)) %>% 
  ggplot(., aes(x=fct_infreq(coding)))+
  geom_bar(stat="count", fill = "dodgerblue4")+
  coord_flip()+
  labs(x="",
       y="Count")+
  theme(panel.spacing = unit(1.5, "lines"),
        text = element_text(size = 12, family = "Times"),
        # axis.text.y = element_text(size=14),
        legend.position = "none",
        panel.grid.major = element_blank(),
        panel.background = element_blank(),
        legend.key=element_blank())

##Table A5: Descriptive Statistics------
disc_table <- full_data %>% 
  mutate(.,
         grade_4 = ifelse(grade == "D",1,0),
         grade_5 = ifelse(grade == "E",1,0),
         grade_6 = ifelse(grade == "F",1,0)) %>% 
  dplyr::select(.,
                male, age, grade_4,grade_5, grade_6)

stargazer(as.data.frame(disc_table),
          covariate.labels = c("Boy", "Age", "Grade 4",
                               "Grade 5", "Grade 6"),
          title = "Descriptive Statistics - Study II",
          label = "tab:dsc_experiment2",
          style = "qje",
          #  notes = c("Therm refers to feelings thermometer, SD refers to social disctance scale, Manip refers to manipulation check."),
          notes.append = T,
          notes.align = "l",
          omit.summary.stat = c("p25", "p75"))



##Table A6: Balance------
balanceTest(treatment~ 
                male + age+
                therm_index_pre +
                therm_arab_pre + therm_imgrnt_pre+ therm_uo_pre +therm_blind_pre +
                cntct_indx_pre  +
                cntct_arab_index_pre + cntct_imgrnt_index_pre + cntct_uo_index_pre + cntct_blind_index_pre+
                similarity_index_pre +
                smlr_arb_pre + smlr_uo_pre + smlr_blind_pre + smlr_imgrnt_pre+
                het_index_pre +
                htr_arab_pre + htr_imgrnt_pre + htr_blind_pre + htr_uo_pre+
                pt_index_pre +
                pt_arab_pre + pt_imgrnt_pre + pt_uo_pre + pt_blind_pre +
                div_scale_pre + cluster(class_id),
              data = full_data)


##Figure A16: Correlations of non-response in endline----
####Thermometor Attrition------
therm_at <- 
  full_data %>% 
  lm_lin(na_therm_index_post ~ 
           treatment,
         covariates = ~  therm_index_pre+
           block_id +male + 
           na_therm_arab_pre + na_therm_blind_pre + 
           na_therm_imgrnt_pre + na_therm_uo_pre,
         clusters = class_id,
         se_type = "stata",
         weights = assignment_prob_weights,
         data = .) 
therm_at_out <-
  tidy(therm_at) %>% 
  filter(.,
         term %in% c("treatment", "therm_index_pre_c")) %>% 
  mutate(.,
         Outcome = "Thermometers \nPost 1",
         term = case_when(
           term == "treatment" ~ "Treatment",
           term == "therm_index_pre_c" ~ "Pre-Treatment \nOutcome",
         ))


therm_at2 <- 
  full_data %>% 
  lm_lin(na_therm_index_post2 ~ 
           treatment,
         covariates = ~  therm_index_pre+
           block_id +male + 
           na_therm_arab_pre + na_therm_blind_pre + 
           na_therm_imgrnt_pre + na_therm_uo_pre,
         clusters = class_id,
         se_type = "stata",
         weights = assignment_prob_weights,
         data = .) 
therm_at_out2 <-
  tidy(therm_at2) %>% 
  filter(.,
         term %in% c("treatment", "therm_index_pre_c")) %>% 
  mutate(.,
         Outcome = "Thermometers \nPost 2",
         term = case_when(
           term == "treatment" ~ "Treatment",
           term == "therm_index_pre_c" ~ "Pre-Treatment \nOutcome",
         ))



###Contact Attrition------

cntct_at <- 
  full_data %>% 
  lm_lin(na_cntct_scale_post ~ 
           treatment,
         covariates = ~  cntct_indx_pre+
           block_id +male +
           na_cntct_1_blind_pre + na_cntct_2_blind_pre +
           na_cntct_3_blind_pre + na_cntct_4_blind_pre +
           na_cntct_1_arab_pre + na_cntct_2_arab_pre + 
           na_cntct_3_arab_pre + na_cntct_4_arab_pre +
           na_cntct_1_imgrnt_pre + na_cntct_2_imgrnt_pre +
           na_cntct_3_imgrnt_pre + na_cntct_4_imgrnt_pre +
           na_cntct_1_uo_pre + na_cntct_2_uo_pre +
           na_cntct_3_uo_pre + na_cntct_4_uo_pre,
         clusters = class_id,
         se_type = "stata",
         weights = assignment_prob_weights,
         data = .) 
cntct_at_out <-
  tidy(cntct_at) %>% 
  filter(.,
         term %in% c("treatment", "cntct_indx_pre_c")) %>% 
  mutate(.,
         Outcome = "Contact Intentions \nPost 1",
         term = case_when(
           term == "treatment" ~ "Treatment",
           term == "cntct_indx_pre_c" ~ "Pre-Treatment \nOutcome",
         ))


cntct_at2 <- 
  full_data %>% 
  lm_lin(na_cntct_scale_post2 ~ 
           treatment,
         covariates = ~  cntct_indx_pre+
           block_id +male +
           na_cntct_1_blind_pre + na_cntct_2_blind_pre +
           na_cntct_3_blind_pre + na_cntct_4_blind_pre +
           na_cntct_1_arab_pre + na_cntct_2_arab_pre + 
           na_cntct_3_arab_pre + na_cntct_4_arab_pre +
           na_cntct_1_imgrnt_pre + na_cntct_2_imgrnt_pre +
           na_cntct_3_imgrnt_pre + na_cntct_4_imgrnt_pre +
           na_cntct_1_uo_pre + na_cntct_2_uo_pre +
           na_cntct_3_uo_pre + na_cntct_4_uo_pre,
         clusters = class_id,
         se_type = "stata",
         weights = assignment_prob_weights,
         data = .) 
cntct_at_out2 <-
  tidy(cntct_at2) %>% 
  filter(.,
         term %in% c("treatment", "cntct_indx_pre_c")) %>% 
  mutate(.,
         Outcome = "Contact Intentions \nPost 2",
         term = case_when(
           term == "treatment" ~ "Treatment",
           term == "cntct_indx_pre_c" ~ "Pre-Treatment \nOutcome",
         ))




###Diveristy Attrition------

div_at <- 
  full_data %>% 
  lm_lin(na_div_scale_post ~ 
           treatment,
         covariates = ~  div_scale_pre+
           block_id +male + 
           na_div1_pre + na_div2_pre + na_div3_pre +
           na_div4_pre + na_div5_pre,
         clusters = class_id,
         se_type = "stata",
         weights = assignment_prob_weights,
         data = .) 
div_at_out <-
  tidy(div_at) %>% 
  filter(.,
         term %in% c("treatment", "div_scale_pre_c")) %>% 
  mutate(.,
         Outcome = "Support for \nDiversity \nPost 1",
         term = case_when(
           term == "treatment" ~ "Treatment",
           term == "div_scale_pre_c" ~ "Pre-Treatment \nOutcome",
         ))


div_at2 <- 
  full_data %>% 
  lm_lin(na_div_scale_post2 ~ 
           treatment,
         covariates = ~  div_scale_pre+
           block_id +male + 
           na_div1_pre + na_div2_pre + na_div3_pre +
           na_div4_pre + na_div5_pre,
         clusters = class_id,
         se_type = "stata",
         weights = assignment_prob_weights,
         data = .) 
div_at_out2 <-
  tidy(div_at2) %>% 
  filter(.,
         term %in% c("treatment", "div_scale_pre_c")) %>% 
  mutate(.,
         Outcome = "Support for \nDiversity \nPost 2",
         term = case_when(
           term == "treatment" ~ "Treatment",
           term == "div_scale_pre_c" ~ "Pre-Treatment \nOutcome",
         ))




###Perspective Taking Attrition------

pt_at <- 
  full_data %>% 
  lm_lin(na_pt_post ~ 
           treatment,
         covariates = ~  pt_index_pre+
           block_id +male +
           na_pt_arab_pre + na_pt_blind_pre +
           na_pt_imgrnt_pre + na_pt_uo_pre,
         clusters = class_id,
         se_type = "stata",
         weights = assignment_prob_weights,
         data = .) 
pt_at_out <-
  tidy(pt_at) %>% 
  filter(.,
         term %in% c("treatment", "pt_index_pre_c")) %>% 
  mutate(.,
         Outcome = "Perspective Taking \nPost 1",
         term = case_when(
           term == "treatment" ~ "Treatment",
           term == "pt_index_pre_c" ~ "Pre-Treatment \nOutcome",
         ))


pt_at2 <- 
  full_data %>% 
  lm_lin(na_pt_post2 ~ 
           treatment,
         covariates = ~  pt_index_pre+
           block_id +male +
           na_pt_arab_pre + na_pt_blind_pre +
           na_pt_imgrnt_pre + na_pt_uo_pre,
         clusters = class_id,
         se_type = "stata",
         weights = assignment_prob_weights,
         data = .) 
pt_at_out2 <-
  tidy(pt_at2) %>% 
  filter(.,
         term %in% c("treatment", "pt_index_pre_c")) %>% 
  mutate(.,
         Outcome = "Perspective Taking \nPost 2",
         term = case_when(
           term == "treatment" ~ "Treatment",
           term == "pt_index_pre_c" ~ "Pre-Treatment \nOutcome",
         ))




###Heterogeneity Attrition------

ht_at <- 
  full_data %>% 
  lm_lin(na_het_post ~ 
           treatment,
         covariates = ~  het_index_pre+
           block_id +male +
           na_htr_arab_pre + na_htr_blind_pre +
           na_htr_imgrnt_pre + na_htr_uo_pre,
         clusters = class_id,
         se_type = "stata",
         weights = assignment_prob_weights,
         data = .) 
ht_at_out <-
  tidy(ht_at) %>% 
  filter(.,
         term %in% c("treatment", "het_index_pre_c")) %>% 
  mutate(.,
         Outcome = "Heterogeneity \nPost 1",
         term = case_when(
           term == "treatment" ~ "Treatment",
           term == "het_index_pre_c" ~ "Pre-Treatment \nOutcome",
         ))

ht_at2 <- 
  full_data %>% 
  lm_lin(na_het_post2 ~ 
           treatment,
         covariates = ~  het_index_pre+
           block_id +male +
           na_htr_arab_pre + na_htr_blind_pre +
           na_htr_imgrnt_pre + na_htr_uo_pre,
         clusters = class_id,
         se_type = "stata",
         weights = assignment_prob_weights,
         data = .) 
ht_at_out2 <-
  tidy(ht_at2) %>% 
  filter(.,
         term %in% c("treatment", "het_index_pre_c")) %>% 
  mutate(.,
         Outcome = "Heterogeneity \nPost 2",
         term = case_when(
           term == "treatment" ~ "Treatment",
           term == "het_index_pre_c" ~ "Pre-Treatment \nOutcome",
         ))


###Similarity Attrition------
sim_at <- 
  full_data %>% 
  lm_lin(na_sim_post ~ 
           treatment,
         covariates = ~  similarity_index_pre+
           block_id +male +
           na_smlr_blind_pre + na_smlr_imgrnt_pre +
           na_smlr_arb_pre + na_smlr_uo_pre,
         clusters = class_id,
         se_type = "stata",
         weights = assignment_prob_weights,
         data = .) 
sim_at_out <-
  tidy(sim_at) %>% 
  filter(.,
         term %in% c("treatment", "similarity_index_pre_c")) %>% 
  mutate(.,
         Outcome = "Similarity \nPost 1",
         term = case_when(
           term == "treatment" ~ "Treatment",
           term == "similarity_index_pre_c" ~ "Pre-Treatment \nOutcome",
         ))



sim_at2 <- 
  full_data %>% 
  lm_lin(na_sim_post2 ~ 
           treatment,
         covariates = ~  similarity_index_pre+
           block_id +male +
           na_smlr_blind_pre + na_smlr_imgrnt_pre +
           na_smlr_arb_pre + na_smlr_uo_pre,
         clusters = class_id,
         se_type = "stata",
         weights = assignment_prob_weights,
         data = .) 
sim_at_out2 <-
  tidy(sim_at2) %>% 
  filter(.,
         term %in% c("treatment", "similarity_index_pre_c")) %>% 
  mutate(.,
         Outcome = "Similarity \nPost 2",
         term = case_when(
           term == "treatment" ~ "Treatment",
           term == "similarity_index_pre_c" ~ "Pre-Treatment \nOutcome",
         ))


####Plot Attrition------

att_models <- 
  rbind(therm_at_out, therm_at_out2, cntct_at_out, cntct_at_out2,
        div_at_out, div_at_out2, pt_at_out, pt_at_out2, ht_at_out, ht_at_out2,
        sim_at_out, sim_at_out2) %>%
  mutate(.,
         se = paste0("(", round(std.error,2),")"),
         pese = paste(round(estimate,2), se),
         term = case_when(
           term == "Treatment" ~ "Correlation of \nTreatment with \nNonresponse",
           term == "Pre-Treatment \nOutcome" ~ "Correlation of \nPre-Treatment \nOutcome with Nonresponse"
         ))
att_models$term <- 
  factor(att_models$term,
         c("Correlation of \nTreatment with \nNonresponse", "Correlation of \nPre-Treatment \nOutcome with Nonresponse"))
att_order<-
  c("Thermometers \nPost 1", "Thermometers \nPost 2", 
    "Contact Intentions \nPost 1", "Contact Intentions \nPost 2",
    "Support for \nDiversity \nPost 1", "Support for \nDiversity \nPost 2",
    "Perspective Taking \nPost 1", "Perspective Taking \nPost 2",
    "Heterogeneity \nPost 1", "Heterogeneity \nPost 2",
    "Similarity \nPost 1", "Similarity \nPost 2") %>% rev()
ggplot(att_models, aes(x=estimate, y=Outcome, color = "dodgerblue4")) + 
  geom_point(position = position_dodge(width = 0.6), show.legend = FALSE)+
  geom_errorbar(aes(xmin=conf.low, xmax=conf.high), 
                width=.1,
                position = position_dodge(width = 0.6),
                show.legend = FALSE) +
  geom_vline(xintercept = 0,
             color = "black",
             linetype = "dashed",
             size = 0.4) +
  xlim(-.3,.3)+
  scale_y_discrete(limits = att_order)+
  labs(x = "",
       y = "",
       color = "")+
  scale_color_manual(values = c("dodgerblue4"))+
  facet_grid(cols = vars(term))+
  # theme_tufte()+
  #guides(color = guide_legend(reverse = TRUE))+
  theme(panel.spacing = unit(1.5, "lines"),
        text = element_text(size = 12, family = "Times"),
        # axis.text.y = element_text(size=14),
        legend.position = "bottom",
        panel.grid.major = element_blank(),
        panel.background = element_blank(),
        legend.key=element_blank())



## Figure A17: ATEs on disagregte measures of contact and affect-----
### Thermometors------

# Arabs
therm_arab <- 
 full_data %>% 
 lm_lin(standardize(therm_arab_post) ~ 
         treatment,
        covariates = ~  therm_arab_pre+
         block_id +male + 
         na_therm_arab_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
therm_arab_out <-
 tidy(therm_arab) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Target = "Arabs",
        Outcome = "Thermometers",
        Wave = "1-2 Weeks")


therm_arab2 <- 
 full_data %>% 
 lm_lin(standardize(therm_arab_post2) ~ 
         treatment,
        covariates = ~  therm_arab_pre+
         block_id +male +
         na_therm_arab_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
therm_arab_out2 <-
 tidy(therm_arab2) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Target = "Arabs",
        Outcome = "Thermometers",
        Wave = "8-13 Weeks")

# Immigrant
therm_imgrnt <- 
 full_data %>% 
 lm_lin(standardize(therm_imgrnt_post) ~ 
         treatment,
        covariates = ~  therm_imgrnt_pre+
         block_id +male + 
         na_therm_imgrnt_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
therm_imgrnt_out <-
 tidy(therm_imgrnt) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Target = "Immigrants",
        Outcome = "Thermometers",
        Wave = "1-2 Weeks")


therm_imgrnt2 <- 
 full_data %>% 
 lm_lin(standardize(therm_imgrnt_post2) ~ 
         treatment,
        covariates = ~  therm_imgrnt_pre+
         block_id +male +
         na_therm_imgrnt_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
therm_imgrnt_out2 <-
 tidy(therm_imgrnt2) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Target = "Immigrants",
        Outcome = "Thermometers",
        Wave = "8-13 Weeks")

# Blind
therm_blind <- 
 full_data %>% 
 lm_lin(standardize(therm_blind_post) ~ 
         treatment,
        covariates = ~  therm_blind_pre+
         block_id +male +
         na_therm_blind_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
therm_blind_out <-
 tidy(therm_blind) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Target = "Visually \nImpaired",
        Outcome = "Thermometers",
        Wave = "1-2 Weeks")

therm_blind2 <- 
 full_data %>% 
 lm_lin(standardize(therm_blind_post2) ~ 
         treatment,
        covariates = ~  therm_blind_pre+
         block_id +male +
         na_therm_blind_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
therm_blind_out2 <-
 tidy(therm_blind2) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Target = "Visually \nImpaired",
        Outcome = "Thermometers",
        Wave = "8-13 Weeks")

therm_uo <- 
 full_data %>% 
 lm_lin(standardize(therm_uo_post) ~ 
         treatment,
        covariates = ~  therm_uo_pre+
         block_id +male +
         na_therm_uo_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
therm_uo_out <-
 tidy(therm_uo) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Target = "Ultra \nOrthodox",
        Outcome = "Thermometers",
        Wave = "1-2 Weeks")


therm_uo2 <- 
 full_data %>% 
 lm_lin(standardize(therm_uo_post2) ~ 
         treatment,
        covariates = ~  therm_uo_pre+
         block_id +male +
         na_therm_uo_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
therm_uo_out2 <-
 tidy(therm_uo2) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Target = "Ultra \nOrthodox",
        Outcome = "Thermometers",
        Wave = "8-13 Weeks")



### Contact Index---------

cntct_arab <- 
 full_data %>% 
 lm_lin(standardize(cntct_arab_index_post) ~ 
         treatment,
        covariates = ~  cntct_arab_index_pre+
         block_id +male + 
         na_cntct_1_arab_pre + na_cntct_2_arab_pre + 
         na_cntct_3_arab_pre + na_cntct_4_arab_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
cntct_arab_out <-
 tidy(cntct_arab) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Outcome = "Contact Intentions",
        Target = "Arabs", 
        Wave = "1-2 Weeks")



cntct_arab2 <- 
 full_data %>% 
 lm_lin(standardize(cntct_arab_index_post2) ~ 
         treatment,
        covariates = ~  cntct_arab_index_pre+
         block_id +male + 
         na_cntct_1_arab_pre + na_cntct_2_arab_pre + 
         na_cntct_3_arab_pre + na_cntct_4_arab_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
cntct_arab_out2 <-
 tidy(cntct_arab2) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Outcome = "Contact Intentions",
        Target = "Arabs", 
        Wave = "8-13 Weeks")


cntct_imgrnt <- 
 full_data %>% 
 lm_lin(standardize(cntct_imgrnt_index_post) ~ 
         treatment,
        covariates = ~  cntct_imgrnt_index_pre+
         block_id +male + 
         na_cntct_1_imgrnt_pre + na_cntct_2_imgrnt_pre +
         na_cntct_3_imgrnt_pre + na_cntct_4_imgrnt_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
cntct_imgrnt_out <-
 tidy(cntct_imgrnt) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Outcome = "Contact Intentions",
        Target = "Immigrants", 
        Wave = "1-2 Weeks")


cntct_imgrnt2 <- 
 full_data %>% 
 lm_lin(standardize(cntct_imgrnt_index_post2) ~ 
         treatment,
        covariates = ~  cntct_imgrnt_index_pre+
         block_id +male+
         na_cntct_1_imgrnt_pre + na_cntct_2_imgrnt_pre +
         na_cntct_3_imgrnt_pre + na_cntct_4_imgrnt_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
cntct_imgrnt_out2 <-
 tidy(cntct_imgrnt2) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Outcome = "Contact Intentions",
        Target = "Immigrants", 
        Wave = "8-13 Weeks")



cntct_blind <- 
 full_data %>% 
 lm_lin(standardize(cntct_blind_index_post) ~ 
         treatment,
        covariates = ~  cntct_blind_index_pre+
         block_id +male +
         na_cntct_1_blind_pre + na_cntct_2_blind_pre +
         na_cntct_3_blind_pre + na_cntct_4_blind_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
cntct_blind_out <-
 tidy(cntct_blind) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Outcome = "Contact Intentions",
        Target = "Visually \nImpaired", 
        Wave = "1-2 Weeks")



cntct_blind2 <- 
 full_data %>% 
 lm_lin(standardize(cntct_blind_index_post2) ~ 
         treatment,
        covariates = ~  cntct_blind_index_pre+
         block_id +male +
         na_cntct_1_blind_pre + na_cntct_2_blind_pre +
         na_cntct_3_blind_pre + na_cntct_4_blind_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
cntct_blind_out2 <-
 tidy(cntct_blind2) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Outcome = "Contact Intentions",
        Target = "Visually \nImpaired", 
        Wave = "8-13 Weeks")




cntct_uo <- 
 full_data %>% 
 lm_lin(standardize(cntct_uo_index_post) ~ 
         treatment,
        covariates = ~  cntct_uo_index_pre+
         block_id +male +
         na_cntct_1_uo_pre + na_cntct_2_uo_pre +
         na_cntct_3_uo_pre + na_cntct_4_uo_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
cntct_uo_out <-
 tidy(cntct_uo) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Outcome = "Contact Intentions",
        Target = "Ultra \nOrthodox", 
        Wave = "1-2 Weeks")



cntct_uo2 <- 
 full_data %>% 
 lm_lin(standardize(cntct_uo_index_post2) ~ 
         treatment,
        covariates = ~  cntct_uo_index_pre+
         block_id +male+
         na_cntct_1_uo_pre + na_cntct_2_uo_pre +
         na_cntct_3_uo_pre + na_cntct_4_uo_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
cntct_uo_out2 <-
 tidy(cntct_uo2) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Outcome = "Contact Intentions",
        Target = "Ultra \nOrthodox", 
        Wave = "8-13 Weeks")





##### Plot Thermometor and Contact disagregated results--------
diss_therm_cntct <-
 rbind(therm_arab_out, therm_arab_out2,
       therm_imgrnt_out, therm_imgrnt_out2,
       therm_blind_out, therm_blind_out2,
       therm_uo_out, therm_uo_out2,
       cntct_arab_out, cntct_arab_out2,
       cntct_imgrnt_out, cntct_imgrnt_out2,
       cntct_blind_out, cntct_blind_out2,
       cntct_uo_out, cntct_uo_out2) %>% 
mutate(.,
       se = paste0("(", round(std.error,2),")"),
       pese = paste(round(estimate,2), se))

diss_therm_cntct$Outcome <- factor(diss_therm_cntct$Outcome,
                       c("Contact Intentions", "Thermometers")) 

diss_therm_cntct$Target <- factor(diss_therm_cntct$Target,
                                   c("Visually \nImpaired","Ultra \nOrthodox",
                                    "Immigrants", "Arabs")) 



ggplot(diss_therm_cntct, aes(x=estimate, y=Target,
                 colour = forcats::fct_rev(Wave))) + 
 geom_point(position = position_dodge(width = 0.6))+
 geom_errorbar(aes(xmin=conf.low, xmax=conf.high), 
               width=.1,
               position = position_dodge(width = 0.6)) +
 geom_vline(xintercept = 0,
            color = "black",
            linetype = "dashed",
            size = 0.4) +
 xlim(-.1,.56)+
 facet_grid(cols = vars(Outcome))+
 geom_text(aes(label= pese, color = forcats::fct_rev(Wave)), 
           position = position_dodge(width = .6),
           hjust = -.82, 
           size = 3,
           family = "Times",
           show.legend = FALSE) +
 labs(x = "",
      y = "",
      color = "")+
 scale_color_manual(values = c("firebrick3","dodgerblue4"))+
 # theme_tufte()+
 guides(color = guide_legend(reverse = TRUE))+
 theme(panel.spacing = unit(1.5, "lines"),
       text = element_text(size = 14, family = "Times"),
       # axis.text.y = element_text(size=14),
       legend.position = "bottom",
       panel.grid.major = element_blank(),
       panel.background = element_blank(),
       legend.key=element_blank())

##Figure A18: ATEs on disaggregated support for diversity-----

div1 <- 
 full_data %>% 
 lm_lin(standardize(div1_post) ~ 
         treatment,
        covariates = ~  div1_pre+
         block_id +male +
         na_div1_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
div1_out <-
 tidy(div1) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Outcome = "Diversity",
        Target = "Important to \nLearn from \nOthers",
        Wave = "1-2 Weeks")



div12 <- 
 full_data %>% 
 lm_lin(standardize(div1_post2) ~ 
         treatment,
        covariates = ~  div1_pre+
         block_id +male + 
         na_div1_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
div1_out2 <-
 tidy(div12) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Outcome = "Diversity",
        Target = "Important to \nLearn from \nOthers",
        Wave = "8-13 Weeks")


div2 <- 
 full_data %>% 
 lm_lin(standardize(div2_post) ~ 
         treatment,
        covariates = ~  div2_pre+
         block_id +male+
         na_div2_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
div2_out <-
 tidy(div2) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Outcome = "Diversity",
        Target = "Important to \nListen to \nOthers",
        Wave = "1-2 Weeks")


div22 <- 
 full_data %>% 
 lm_lin(standardize(div2_post2) ~ 
         treatment,
        covariates = ~  div2_pre+
         block_id +male+
         na_div2_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
div2_out2 <-
 tidy(div22) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Outcome = "Diversity",
        Target = "Important to \nListen to \nOthers",
        Wave = "8-13 Weeks")



div3 <- 
 full_data %>% 
 lm_lin(standardize(div3_post) ~ 
         treatment,
        covariates = ~  div3_pre+
         block_id +male+
         na_div3_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
div3_out <-
 tidy(div3) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Outcome = "Diversity",
        Target = "Enjoy Learning \nWith Others",
        Wave = "1-2 Weeks")



div32 <- 
 full_data %>% 
 lm_lin(standardize(div3_post2) ~ 
         treatment,
        covariates = ~  div3_pre+
         block_id +male+
         na_div3_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
div3_out2 <-
 tidy(div32) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Outcome = "Diversity",
        Target = "Enjoy Learning \nWith Others",
        Wave = "8-13 Weeks")


div4 <- 
 full_data %>% 
 lm_lin(standardize(div4_post) ~ 
         treatment,
        covariates = ~  div4_pre+
         block_id +male+
         na_div4_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
div4_out <-
 tidy(div4) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Outcome = "Diversity",
        Target = "Enjoy Playing \nWith Others",
        Wave = "1-2 Weeks")



div42 <- 
 full_data %>% 
 lm_lin(standardize(div4_post2) ~ 
         treatment,
        covariates = ~  div4_pre+
         block_id +male+
         na_div4_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
div4_out2 <-
 tidy(div42) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Outcome = "Diversity",
        Target = "Enjoy Playing \nWith Others",
        Wave = "8-13 Weeks")


div5<- 
 full_data %>% 
 lm_lin(standardize(div5_post) ~ 
         treatment,
        covariates = ~  div5_pre+
         block_id +male+
         na_div5_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
div5_out <-
 tidy(div5) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Outcome = "Diversity",
        Target = "Diversity Helps \nTeamwork",
        Wave = "1-2 Weeks")



div52<- 
 full_data %>% 
 lm_lin(standardize(div5_post2) ~ 
         treatment,
        covariates = ~  div5_pre+
         block_id +male+
         na_div5_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
div5_out2 <-
 tidy(div52) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Outcome = "Diversity",
        Target = "Diversity Helps \nTeamwork",
        Wave = "8-13 Weeks")



#### Plot Diversity effects disaggerageted------


div_disag <- rbind(div1_out, div1_out2,
                   div2_out, div2_out2,
                   div3_out, div3_out2,
                   div4_out, div4_out2,
                   div5_out, div5_out2) %>% 
 mutate(.,
        se = paste0("(", round(std.error,2),")"),
        pese = paste(round(estimate,2), se))


annotation_mech <- data.frame(
 x = c(-0.1,-0.1),
 y = c(4.15,3.85),
 label = c("1-2 Weeks", "8-13 Weeks")
)

ggplot(div_disag, aes(x=estimate, y=Target,
                 colour = forcats::fct_rev(Wave))) + 
 geom_point(position = position_dodge(width = 0.6))+
 geom_errorbar(aes(xmin=conf.low, xmax=conf.high), 
               width=.1,
               position = position_dodge(width = 0.6)) +
 geom_vline(xintercept = 0,
            color = "black",
            linetype = "dashed",
            size = 0.4) +
 xlim(-.18,.57)+
 geom_text(aes(label= pese, color = forcats::fct_rev(Wave)), 
           position = position_dodge(width = .6),
           hjust = -1.8, 
           size = 3,
           family = "Times",
           show.legend = FALSE) +
 geom_label(data=annotation_mech, aes(x=x, y=y, label=label),
            color=c("dodgerblue4", "firebrick3"), 
            size=2 , angle=45, fontface="bold")+
 labs(x = "",
      y = "",
      color = "")+
 scale_color_manual(values = c("firebrick3","dodgerblue4"))+
 # theme_tufte()+
 guides(color = guide_legend(reverse = TRUE))+
 theme(panel.spacing = unit(1.5, "lines"),
       text = element_text(size = 14, family = "Times"),
       # axis.text.y = element_text(size=14),
       legend.position = "none",
       panel.grid.major = element_blank(),
       panel.background = element_blank(),
       legend.key=element_blank())



##Figure A19: ATEs on group-specific mechanism measures------
###Perspective Taking-----

pt_arabs <- 
 full_data %>% 
 lm_lin(standardize(pt_arab_post) ~ 
         treatment,
        covariates = ~  pt_arab_pre+
         block_id +male + na_pt_arab_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
pt_arabs_out <-
 tidy(pt_arabs) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Target = "Arabs",
        Outcome = "Perspective Taking",
        Wave = "1-2 Weeks")


pt_arabs2 <- 
 full_data %>% 
 lm_lin(standardize(pt_arab_post2) ~ 
         treatment,
        covariates = ~  pt_arab_pre+
         block_id +male + na_pt_arab_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
pt_arabs_out2 <-
 tidy(pt_arabs2) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Target = "Arabs",
        Outcome = "Perspective Taking",
        Wave = "8-13 Weeks")


pt_imgrnt <- 
 full_data %>% 
 lm_lin(standardize(pt_imgrnt_post) ~ 
         treatment,
        covariates = ~  pt_imgrnt_pre+
         block_id +male + na_pt_imgrnt_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
pt_imgrnt_out <-
 tidy(pt_imgrnt) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Target = "Immigrants",
        Outcome = "Perspective Taking",
        Wave = "1-2 Weeks")



pt_imgrnt2 <- 
 full_data %>% 
 lm_lin(standardize(pt_imgrnt_post2) ~ 
         treatment,
        covariates = ~  pt_imgrnt_pre+
         block_id +male + na_pt_imgrnt_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
pt_imgrnt_out2 <-
 tidy(pt_imgrnt2) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Target = "Immigrants",
        Outcome = "Perspective Taking",
        Wave = "8-13 Weeks")


pt_blind <- 
 full_data %>% 
 lm_lin(standardize(pt_blind_post) ~ 
         treatment,
        covariates = ~  pt_blind_pre+
         block_id +male + na_pt_blind_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
pt_blind_out <-
 tidy(pt_blind) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Target = "Visually \nImpaired",
        Outcome = "Perspective Taking",
        Wave = "1-2 Weeks")



pt_blind2 <- 
 full_data %>% 
 lm_lin(standardize(pt_blind_post2) ~ 
         treatment,
        covariates = ~  pt_blind_pre+
         block_id +male + na_pt_blind_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
pt_blind_out2 <-
 tidy(pt_blind2) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Target = "Visually \nImpaired",
        Outcome = "Perspective Taking",
        Wave = "8-13 Weeks")



pt_uo <- 
 full_data %>% 
 lm_lin(standardize(pt_uo_post) ~ 
         treatment,
        covariates = ~  pt_uo_pre+
         block_id +male + na_pt_uo_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
pt_uo_out <-
 tidy(pt_uo) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Target = "Ultra \nOrthodox",
        Outcome = "Perspective Taking",
        Wave = "1-2 Weeks")



pt_uo2 <- 
 full_data %>% 
 lm_lin(standardize(pt_uo_post2) ~ 
         treatment,
        covariates = ~  pt_uo_pre+
         block_id +male + na_pt_uo_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
pt_uo_out2 <-
 tidy(pt_uo2) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Target = "Ultra \nOrthodox",
        Outcome = "Perspective Taking",
        Wave = "8-13 Weeks")


###Heterogeneity-------
het_arabs <- 
 full_data %>% 
 lm_lin(standardize(htr_arab_post) ~ 
         treatment,
        covariates = ~  htr_arab_pre+
         block_id +male + na_htr_arab_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
het_arabs_out <-
 tidy(het_arabs) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Target = "Arabs",
        Outcome = "Group Heterogeneity",
        Wave = "1-2 Weeks")


het_arabs2 <- 
 full_data %>% 
 lm_lin(standardize(htr_arab_post2) ~ 
         treatment,
        covariates = ~  htr_arab_pre+
         block_id +male + na_htr_arab_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
het_arabs_out2 <-
 tidy(het_arabs2) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Target = "Arabs",
        Outcome = "Group Heterogeneity",
        Wave = "8-13 Weeks")




het_imgrnt <- 
 full_data %>% 
 lm_lin(standardize(htr_imgrnt_post) ~ 
         treatment,
        covariates = ~  htr_imgrnt_pre+
         block_id +male + na_htr_imgrnt_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
het_imgrnt_out <-
 tidy(het_imgrnt) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Target = "Immigrants",
        Outcome = "Group Heterogeneity",
        Wave = "1-2 Weeks")




het_imgrnt2 <- 
 full_data %>% 
 lm_lin(standardize(htr_imgrnt_post2) ~ 
         treatment,
        covariates = ~  htr_imgrnt_pre+
         block_id +male + na_htr_imgrnt_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
het_imgrnt_out2 <-
 tidy(het_imgrnt2) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Target = "Immigrants",
        Outcome = "Group Heterogeneity",
        Wave = "8-13 Weeks")


het_blind <- 
 full_data %>% 
 lm_lin(standardize(htr_blind_post) ~ 
         treatment,
        covariates = ~  htr_blind_pre+
         block_id +male + na_htr_blind_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
het_blind_out <-
 tidy(het_blind) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Target = "Visually \nImpaired",
        Outcome = "Group Heterogeneity",
        Wave = "1-2 Weeks")



het_blind2 <- 
 full_data %>% 
 lm_lin(standardize(htr_blind_post2) ~ 
         treatment,
        covariates = ~  htr_blind_pre+
         block_id +male + na_htr_blind_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
het_blind_out2 <-
 tidy(het_blind2) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Target = "Visually \nImpaired",
        Outcome = "Group Heterogeneity",
        Wave = "8-13 Weeks")



het_uo <- 
 full_data %>% 
 lm_lin(standardize(htr_uo_post) ~ 
         treatment,
        covariates = ~  htr_uo_pre+
         block_id +male + na_htr_uo_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
het_uo_out <-
 tidy(het_uo) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Target = "Ultra \nOrthodox",
        Outcome = "Group Heterogeneity",
        Wave = "1-2 Weeks")



het_uo2 <- 
 full_data %>% 
 lm_lin(standardize(htr_uo_post2) ~ 
         treatment,
        covariates = ~  htr_uo_pre+
         block_id +male + na_htr_uo_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
het_uo_out2 <-
 tidy(het_uo2) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Target = "Ultra \nOrthodox",
        Outcome = "Group Heterogeneity",
        Wave = "8-13 Weeks")


### Similarity----
sim_arabs <- 
 full_data %>% 
 lm_lin(standardize(smlr_arb_post) ~ 
         treatment,
        covariates = ~  smlr_arb_pre+
         block_id +male + na_smlr_arb_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
sim_arabs_out <-
 tidy(sim_arabs) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Target = "Arabs",
        Outcome = "Group Similarity",
        Wave = "1-2 Weeks")


sim_arabs2 <- 
 full_data %>% 
 lm_lin(standardize(smlr_arb_post2) ~ 
         treatment,
        covariates = ~  smlr_arb_pre+
         block_id +male + na_smlr_arb_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
sim_arabs_out2 <-
 tidy(sim_arabs2) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Target = "Arabs",
        Outcome = "Group Similarity",
        Wave = "8-13 Weeks")




sim_imgrnt <- 
 full_data %>% 
 lm_lin(standardize(smlr_imgrnt_post) ~ 
         treatment,
        covariates = ~  smlr_imgrnt_pre+
         block_id +male + na_smlr_imgrnt_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
sim_imgrnt_out <-
 tidy(sim_imgrnt) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Target = "Immigrants",
        Outcome = "Group Similarity",
        Wave = "1-2 Weeks")



sim_imgrnt2 <- 
 full_data %>% 
 lm_lin(standardize(smlr_imgrnt_post2) ~ 
         treatment,
        covariates = ~  smlr_imgrnt_pre+
         block_id +male + na_smlr_imgrnt_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
sim_imgrnt_out2 <-
 tidy(sim_imgrnt2) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Target = "Immigrants",
        Outcome = "Group Similarity",
        Wave = "8-13 Weeks")


sim_blind <- 
 full_data %>% 
 lm_lin(standardize(smlr_blind_post) ~ 
         treatment,
        covariates = ~  smlr_blind_pre+
         block_id +male + na_smlr_blind_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
sim_blind_out <-
 tidy(sim_blind) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Target = "Visually \nImpaired",
        Outcome = "Group Similarity",
        Wave = "1-2 Weeks")


sim_blind2 <- 
 full_data %>% 
 lm_lin(standardize(smlr_blind_post2) ~ 
         treatment,
        covariates = ~  smlr_blind_pre+
         block_id +male + na_smlr_blind_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
sim_blind_out2 <-
 tidy(sim_blind2) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Target = "Visually \nImpaired",
        Outcome = "Group Similarity",
        Wave = "8-13 Weeks")



sim_uo <- 
 full_data %>% 
 lm_lin(standardize(smlr_uo_post) ~ 
         treatment,
        covariates = ~  smlr_uo_pre+
         block_id +male + na_smlr_uo_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
sim_uo_out <-
 tidy(sim_uo) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Target = "Ultra \nOrthodox",
        Outcome = "Group Similarity",
        Wave = "1-2 Weeks")



sim_uo2 <- 
 full_data %>% 
 lm_lin(standardize(smlr_uo_post2) ~ 
         treatment,
        covariates = ~  smlr_uo_pre+
         block_id +male + na_smlr_uo_pre,
        clusters = class_id,
        se_type = "stata",
        weights = assignment_prob_weights,
        data = .) 
sim_uo_out2 <-
 tidy(sim_uo2) %>% 
 filter(.,
        term == "treatment") %>% 
 mutate(.,
        Target = "Ultra \nOrthodox",
        Outcome = "Group Similarity",
        Wave = "8-13 Weeks")



##### Plot group specific mechanisms--------
diss_mech <-
 rbind(het_arabs_out, het_arabs_out2, het_imgrnt_out, het_imgrnt_out2, het_uo_out, het_uo_out2, het_blind_out, het_blind_out2,
       sim_arabs_out, sim_arabs_out2, sim_imgrnt_out, sim_imgrnt_out2, sim_uo_out, sim_uo_out2, sim_blind_out, sim_blind_out2,
       pt_arabs_out, pt_arabs_out2, pt_imgrnt_out, pt_imgrnt_out2, pt_uo_out, pt_uo_out2, pt_blind_out, pt_blind_out2) %>% 
 mutate(.,
        se = paste0("(", round(std.error,2),")"),
        pese = paste(round(estimate,2), se))

diss_mech$Outcome <- factor(diss_mech$Outcome,
                                   c("Group Heterogeneity","Group Similarity", "Perspective Taking")) 

diss_mech$Target <- factor(diss_mech$Target,
                                  c("Visually \nImpaired","Ultra \nOrthodox",
                                    "Immigrants", "Arabs")) 



ggplot(diss_mech, aes(x=estimate, y=Target,
                             colour = forcats::fct_rev(Wave))) + 
 geom_point(position = position_dodge(width = 0.6))+
 geom_errorbar(aes(xmin=conf.low, xmax=conf.high), 
               width=.1,
               position = position_dodge(width = 0.6)) +
 geom_vline(xintercept = 0,
            color = "black",
            linetype = "dashed",
            size = 0.4) +
 xlim(-.25,.8)+
 facet_grid(cols = vars(Outcome))+
 geom_text(aes(label= pese, color = forcats::fct_rev(Wave)), 
           position = position_dodge(width = .6),
           hjust = -.8, 
           size = 3,
           family = "Times",
           show.legend = FALSE) +
 labs(x = "",
      y = "",
      color = "")+
 scale_color_manual(values = c("firebrick3","dodgerblue4"))+
 # theme_tufte()+
 guides(color = guide_legend(reverse = TRUE))+
 theme(panel.spacing = unit(1.5, "lines"),
       text = element_text(size = 12, family = "Times"),
       # axis.text.y = element_text(size=14),
       legend.position = "bottom",
       panel.grid.major = element_blank(),
       panel.background = element_blank(),
       legend.key=element_blank())


##Figure A20: Pre-treatment levels of perspective taking-----
# select relevant measures
pt_data <- full_data %>% 
  dplyr::select(.,
                pt_arab_pre, pt_uo_pre,
                pt_blind_pre, pt_imgrnt_pre)

# pivot data
pt_piv <-
  pt_data %>% 
  pivot_longer(pt_arab_pre:pt_imgrnt_pre, 
               names_to = "group", 
               values_to = "value") 

# get mean & SD
pt_estimate <- 
  pt_piv %>%
  group_by(group) %>%
  summarize(Mean_pt = mean(value, na.rm=TRUE),
            Sd_pt = sd(value, na.rm=TRUE),
            n_pt = n()) %>% 
  mutate(se_pt = Sd_pt / sqrt(n_pt),
         lower_ci_rust = Mean_pt - qt(1 - (0.05 / 2), n_pt- 1) * se_pt,
         upper_ci_rust = Mean_pt + qt(1 - (0.05 / 2), n_pt - 1) * se_pt) %>% 
  mutate(.,
         group = case_when(
           group == "pt_blind_pre" ~ "Visually \nImpaired",
           group == "pt_arab_pre" ~ "Arab",
           group == "pt_uo_pre" ~ "Ultra-Orthodox",
           group == "pt_imgrnt_pre" ~ "Immigrant"
         ))

ggplot(pt_estimate, aes(x=Mean_pt, y=group, colour="dodgerblue4")) + 
  geom_errorbar(aes(xmin=lower_ci_rust, xmax=upper_ci_rust), width=.1) +
  geom_point()+
  labs(x = "Important to Take Outgroup Perspective",
       y = "")+
  scale_color_manual(values = c("dodgerblue4"))+
  geom_vline(xintercept = mean(pt_piv$value, na.rm = T),
             color = "black",
             linetype = "dashed",
             size = 0.4)+
  xlim(1,5)+
  theme(text = element_text(size = 14, family = "Times"),
        legend.position = "none",
        panel.grid.major = element_blank(),
        panel.background = element_blank())


##Figure A21: Interaction Diagnostic------
# Diagnose common support and linearity
as.data.frame(full_data) %>% 
  interflex(Y = "therm_index_post", 
            D = "treatment", 
            X = "therm_index_pre", data =., 
            estimator = "raw", 
            na.rm = T,
            Ylabel = "Post-Treatment Thermometer", Dlabel = "Treatment", 
            Xlabel="Pre-Treatment Thermometer", main = "",
            theme.bw = T, show.grid = F,
            cex.main = 1.2, ncols=2)



##Figure A22: Moderation by pre-treatment therm------
as.data.frame(full_data) %>% 
  mutate(.,
         block_id = as.factor(block_id)) %>% 
  interflex(Y = "therm_index_post",
            D = "treatment", X = "therm_index_pre", 
            Z = c("block_id", "male", "na_therm_arab_pre", "na_therm_blind_pre", 
                  "na_therm_imgrnt_pre", "na_therm_uo_pre"),
            nboots = 1000, nsimu = 1000,
            data =., estimator = "binning",  vcov.type = "cluster", 
            nbins = 3, cl = "class_id", na.rm = T, main = "Marginal Effects") %>% 
plot(.,
     theme.bw = T, show.grid = F,
     Ylabel = "Post Thermometer", 
     Dlabel = "Treatment", Xlabel = "Pre-Treatment Thermometer",
     cex.lab = .5)


##Figure A23: ATEs only among respondents to both waves-----
###Thermometor------
o_therm_index <- 
  full_data %>% 
  filter(.,
         !is.na(therm_index_post),
         !is.na(therm_index_post2)) %>% 
  lm_lin(standardize(therm_index_post) ~ 
           treatment,
         covariates = ~  therm_index_pre+
           block_id +male + 
           na_therm_arab_pre + na_therm_blind_pre + 
           na_therm_imgrnt_pre + na_therm_uo_pre,
         clusters = class_id,
         se_type = "stata",
         weights = assignment_prob_weights,
         data = .) 
o_therm_index_out <-
  tidy(o_therm_index) %>% 
  filter(.,
         term == "treatment") %>% 
  mutate(.,
         Target = "Index",
         Outcome = "Thermometers",
         Wave = "1-2 Weeks")


o_therm_index2 <- 
  full_data %>% 
  filter(.,
         !is.na(therm_index_post),
         !is.na(therm_index_post2)) %>% 
  lm_lin(standardize(therm_index_post2) ~ 
           treatment,
         covariates = ~  therm_index_pre+
           block_id +male + 
           na_therm_arab_pre + na_therm_blind_pre + 
           na_therm_imgrnt_pre + na_therm_uo_pre,
         clusters = class_id,
         se_type = "stata",
         weights = assignment_prob_weights,
         data = .) 
o_therm_index_out2 <-
  tidy(o_therm_index2) %>% 
  filter(.,
         term == "treatment") %>% 
  mutate(.,
         Target = "Index",
         Outcome = "Thermometers",
         Wave = "8-13 Weeks")

###Contact------
o_cntct_index <- 
  full_data %>% 
  filter(.,
         !is.na(cntct_indx_post),
         !is.na(cntct_indx_post2)) %>% 
  lm_lin(standardize(cntct_indx_post) ~ 
           treatment,
         covariates = ~  cntct_indx_pre+
           block_id +male +
           na_cntct_1_blind_pre + na_cntct_2_blind_pre +
           na_cntct_3_blind_pre + na_cntct_4_blind_pre +
           na_cntct_1_arab_pre + na_cntct_2_arab_pre + 
           na_cntct_3_arab_pre + na_cntct_4_arab_pre +
           na_cntct_1_imgrnt_pre + na_cntct_2_imgrnt_pre +
           na_cntct_3_imgrnt_pre + na_cntct_4_imgrnt_pre +
           na_cntct_1_uo_pre + na_cntct_2_uo_pre +
           na_cntct_3_uo_pre + na_cntct_4_uo_pre,
         clusters = class_id,
         se_type = "stata",
         weights = assignment_prob_weights,
         data = .) 
o_cntct_index_out <-
  tidy(o_cntct_index) %>% 
  filter(.,
         term == "treatment") %>% 
  mutate(.,
         Outcome = "Contact \nIntentions",
         Target = "Index",
         Wave = "1-2 Weeks")


o_cntct_index2 <- 
  full_data %>% 
  filter(.,
         !is.na(cntct_indx_post),
         !is.na(cntct_indx_post2)) %>% 
  lm_lin(standardize(cntct_indx_post2) ~ 
           treatment,
         covariates = ~  cntct_indx_pre+
           block_id +male+ 
           na_cntct_1_blind_pre + na_cntct_2_blind_pre +
           na_cntct_3_blind_pre + na_cntct_4_blind_pre +
           na_cntct_1_arab_pre + na_cntct_2_arab_pre + 
           na_cntct_3_arab_pre + na_cntct_4_arab_pre +
           na_cntct_1_imgrnt_pre + na_cntct_2_imgrnt_pre +
           na_cntct_3_imgrnt_pre + na_cntct_4_imgrnt_pre +
           na_cntct_1_uo_pre + na_cntct_2_uo_pre +
           na_cntct_3_uo_pre + na_cntct_4_uo_pre,
         clusters = class_id,
         se_type = "stata",
         weights = assignment_prob_weights,
         data = .) 
o_cntct_index_out2 <-
  tidy(o_cntct_index2) %>% 
  filter(.,
         term == "treatment") %>% 
  mutate(.,
         Outcome = "Contact \nIntentions",
         Target = "Index",
         Wave = "8-13 Weeks")

###Contact Behavior------
o_cntct_bhvr <- 
  full_data %>% 
  filter(.,
         !is.na(contact_project),
         !is.na(contact_project_post2)) %>% 
  lm_lin(standardize(contact_project) ~ 
           treatment,
         covariates = ~  cntct_indx_pre+ therm_index_pre + div_scale_pre+
           block_id +male +
           na_cntct_1_blind_pre + na_cntct_2_blind_pre +
           na_cntct_3_blind_pre + na_cntct_4_blind_pre +
           na_cntct_1_arab_pre + na_cntct_2_arab_pre + 
           na_cntct_3_arab_pre + na_cntct_4_arab_pre +
           na_cntct_1_imgrnt_pre + na_cntct_2_imgrnt_pre +
           na_cntct_3_imgrnt_pre + na_cntct_4_imgrnt_pre +
           na_cntct_1_uo_pre + na_cntct_2_uo_pre +
           na_cntct_3_uo_pre + na_cntct_4_uo_pre+
           na_therm_arab_pre + na_therm_blind_pre + 
           na_therm_imgrnt_pre + na_therm_uo_pre+
           na_div1_pre + na_div2_pre + na_div3_pre +
           na_div4_pre + na_div5_pre,
         clusters = class_id,
         se_type = "stata",
         weights = assignment_prob_weights,
         data = .) 
o_cntct_bhvr_out <-
  tidy(o_cntct_bhvr) %>% 
  filter(.,
         term == "treatment") %>% 
  mutate(.,
         Outcome = "Register \nContact",
         Target = "Behavioral \nMeasure",
         Wave = "1-2 Weeks")



o_cntct_bhvr2 <- 
  full_data %>% 
  filter(.,
         !is.na(contact_project),
         !is.na(contact_project_post2)) %>% 
  lm_lin(standardize(contact_project_post2) ~ 
           treatment,
         covariates = ~  cntct_indx_pre+ therm_index_pre + div_scale_pre+ 
           block_id +male + 
           na_cntct_1_blind_pre + na_cntct_2_blind_pre +
           na_cntct_3_blind_pre + na_cntct_4_blind_pre +
           na_cntct_1_arab_pre + na_cntct_2_arab_pre + 
           na_cntct_3_arab_pre + na_cntct_4_arab_pre +
           na_cntct_1_imgrnt_pre + na_cntct_2_imgrnt_pre +
           na_cntct_3_imgrnt_pre + na_cntct_4_imgrnt_pre +
           na_cntct_1_uo_pre + na_cntct_2_uo_pre +
           na_cntct_3_uo_pre + na_cntct_4_uo_pre+
           na_therm_arab_pre + na_therm_blind_pre + 
           na_therm_imgrnt_pre + na_therm_uo_pre+
           na_div1_pre + na_div2_pre + na_div3_pre +
           na_div4_pre + na_div5_pre,
         clusters = class_id,
         se_type = "stata",
         weights = assignment_prob_weights,
         data = .) 
o_cntct_bhvr_out2 <-
  tidy(o_cntct_bhvr2) %>% 
  filter(.,
         term == "treatment") %>% 
  mutate(.,
         Outcome = "Register \nContact",
         Target = "Behavioral \nMeasure",
         Wave = "8-13 Weeks")




###Diversity------
o_div_index <- 
  full_data %>% 
  filter(.,
         !is.na(div_scale_post),
         !is.na(div_scale_post2)) %>% 
  lm_lin(standardize(div_scale_post) ~ 
           treatment,
         covariates = ~  div_scale_pre+
           block_id +male + 
           na_div1_pre + na_div2_pre + na_div3_pre +
           na_div4_pre + na_div5_pre,
         clusters = class_id,
         se_type = "stata",
         weights = assignment_prob_weights,
         data = .) 
o_div_index_out <-
  tidy(o_div_index) %>% 
  filter(.,
         term == "treatment") %>% 
  mutate(.,
         Outcome = "Diversity \nAttitudes",
         Target = "Index",
         Wave = "1-2 Weeks")

o_div_index2 <- 
  full_data %>% 
  filter(.,
         !is.na(div_scale_post),
         !is.na(div_scale_post2)) %>% 
  lm_lin(standardize(div_scale_post2) ~ 
           treatment,
         covariates = ~  div_scale_pre+
           block_id +male + 
           na_div1_pre + na_div2_pre + na_div3_pre +
           na_div4_pre + na_div5_pre,
         clusters = class_id,
         se_type = "stata",
         weights = assignment_prob_weights,
         data = .) 
o_div_index_out2 <-
  tidy(o_div_index2) %>% 
  filter(.,
         term == "treatment") %>% 
  mutate(.,
         Outcome = "Diversity \nAttitudes",
         Target = "Index",
         Wave = "8-13 Weeks")



### Diversity Behavior------

o_div_bhvr2 <- 
  full_data %>% 
  filter(.,
         school != "Democrati") %>% 
  lm_lin(standardize(wrist_band) ~ 
           treatment,
         covariates = ~  cntct_indx_pre+ therm_index_pre + div_scale_pre+ 
           block_id + male +
           na_cntct_1_blind_pre + na_cntct_2_blind_pre +
           na_cntct_3_blind_pre + na_cntct_4_blind_pre +
           na_cntct_1_arab_pre + na_cntct_2_arab_pre + 
           na_cntct_3_arab_pre + na_cntct_4_arab_pre +
           na_cntct_1_imgrnt_pre + na_cntct_2_imgrnt_pre +
           na_cntct_3_imgrnt_pre + na_cntct_4_imgrnt_pre +
           na_cntct_1_uo_pre + na_cntct_2_uo_pre +
           na_cntct_3_uo_pre + na_cntct_4_uo_pre+
           na_therm_arab_pre + na_therm_blind_pre + 
           na_therm_imgrnt_pre + na_therm_uo_pre+
           na_div1_pre + na_div2_pre + na_div3_pre +
           na_div4_pre + na_div5_pre,
         clusters = class_id,
         se_type = "stata",
         weights = assignment_prob_weights,
         data = .) 
o_div_bhvr_out2 <-
  tidy(o_div_bhvr2) %>% 
  filter(.,
         term == "treatment") %>% 
  mutate(.,
         Outcome = "Diversity \nWristband",
         Target = "Behavioral \nMeasure",
         Wave = "8-13 Weeks")


#### Plot results------

summary_study2 <- rbind(o_therm_index_out, o_therm_index_out2,
                        o_cntct_index_out, o_cntct_index_out2,
                        o_cntct_bhvr_out, o_cntct_bhvr_out2,
                        o_div_index_out, o_div_index_out2,
                        o_div_bhvr_out2) %>% 
  mutate(.,
         se = paste0("(", round(std.error,2),")"),
         pese = paste(round(estimate,2), se))

# set order
summary_study2$Outcome <- 
  factor(summary_study2$Outcome, 
         c("Diversity \nWristband", "Register \nContact", 
           "Diversity \nAttitudes",
           "Contact \nIntentions",
           "Thermometers"))

# Add annotation
annotation <- data.frame(
  x = c(-0.1,-0.1),
  y = c(5.2,4.8),
  label = c("1-2 Weeks", "8-13 Weeks")
)

ggplot(summary_study2, aes(x=estimate, y=Outcome,
                           colour = forcats::fct_rev(Wave))) + 
  geom_point(position = position_dodge(width = 0.6))+
  geom_errorbar(aes(xmin=conf.low, xmax=conf.high),
                position = position_dodge(width = 0.6),
                width=.1) +
  geom_vline(xintercept = 0,
             color = "black",
             linetype = "dashed",
             size = 0.4) +
  xlim(-0.18,.7)+
  geom_text(aes(label= pese, color = forcats::fct_rev(Wave)), 
            position = position_dodge(width = .6),
            family = "Times",
            hjust = -1, 
            size = 3,
            show.legend = FALSE) +
  geom_label(data=annotation, aes(x=x, y=y, label=label),
             color=c("dodgerblue4", "firebrick3"), 
             size=2 , angle=45, fontface="bold")+
  labs(x = "",
       y = "",
       color = "")+
  scale_color_manual(values = c("firebrick3", "dodgerblue4"))+
  guides(color = guide_legend(reverse = TRUE))+
  theme(panel.spacing = unit(1.5, "lines"),
        text = element_text(size = 15, family = "Times"),
        legend.position = "none",
        panel.grid.major = element_blank(),
        panel.background = element_blank(),
        legend.key=element_blank())



##Figure A24:Mechanisms among respondents to both waves------

###Perspective Taking-----
o_pt_index <- 
  full_data %>% 
  filter(.,
         !is.na(pt_index_post)) %>% 
  lm_lin(standardize(pt_index_post) ~ 
           treatment,
         covariates = ~  pt_index_pre+
           block_id +male +
           na_pt_arab_pre + na_pt_blind_pre +
           na_pt_imgrnt_pre + na_pt_uo_pre,
         clusters = class_id,
         se_type = "stata",
         weights = assignment_prob_weights,
         data = .) 
o_pt_index_out <-
  tidy(o_pt_index) %>% 
  filter(.,
         term == "treatment") %>% 
  mutate(.,
         Target = "Index",
         Outcome = "Perspective Taking",
         Wave = "1-2 Weeks")


o_pt_index2 <- 
  full_data %>% 
  filter(.,
         !is.na(pt_index_post2)) %>% 
  lm_lin(standardize(pt_index_post2) ~ 
           treatment,
         covariates = ~  pt_index_pre+
           block_id +male +
           na_pt_arab_pre + na_pt_blind_pre +
           na_pt_imgrnt_pre + na_pt_uo_pre,
         clusters = class_id,
         se_type = "stata",
         weights = assignment_prob_weights,
         data = .) 
o_pt_index_out2 <-
  tidy(o_pt_index2) %>% 
  filter(.,
         term == "treatment") %>% 
  mutate(.,
         Target = "Index",
         Outcome = "Perspective Taking",
         Wave = "8-13 Weeks")

### Heterogeneity------

o_het_index <- 
  full_data %>% 
  filter(.,
         !is.na(het_index_post)) %>% 
  lm_lin(standardize(het_index_post) ~ 
           treatment,
         covariates = ~  het_index_pre+
           block_id +male +
           na_htr_arab_pre + na_htr_blind_pre +
           na_htr_imgrnt_pre + na_htr_uo_pre,
         clusters = class_id,
         se_type = "stata",
         weights = assignment_prob_weights,
         data = .) 
o_het_index_out <-
  tidy(o_het_index) %>% 
  filter(.,
         term == "treatment") %>% 
  mutate(.,
         Target = "Index",
         Outcome = "Group Heterogeneity",
         Wave = "1-2 Weeks")



o_het_index2 <- 
  full_data %>% 
  filter(.,
         !is.na(het_index_post2)) %>% 
  lm_lin(standardize(het_index_post2) ~ 
           treatment,
         covariates = ~  het_index_pre+
           block_id +male +
           na_htr_arab_pre + na_htr_blind_pre +
           na_htr_imgrnt_pre + na_htr_uo_pre,
         clusters = class_id,
         se_type = "stata",
         weights = assignment_prob_weights,
         data = .) 
o_het_index_out2 <-
  tidy(o_het_index2) %>% 
  filter(.,
         term == "treatment") %>% 
  mutate(.,
         Target = "Index",
         Outcome = "Group Heterogeneity",
         Wave = "8-13 Weeks")



### Similarity------

o_sim_index <- 
  full_data %>%
  filter(.,
         !is.na(similarity_index_post)) %>% 
  filter(.,
         !is.na(similarity_index_post)) %>% 
  lm_lin(standardize(similarity_index_post) ~ 
           treatment,
         covariates = ~  similarity_index_pre+
           block_id +male +
           na_smlr_blind_pre + na_smlr_imgrnt_pre +
           na_smlr_arb_pre + na_smlr_uo_pre,
         clusters = class_id,
         se_type = "stata",
         weights = assignment_prob_weights,
         data = .) 
o_sim_index_out <-
  tidy(o_sim_index) %>% 
  filter(.,
         term == "treatment") %>% 
  mutate(.,
         Target = "Index",
         Outcome = "Group Similarity",
         Wave = "1-2 Weeks")


o_sim_index2 <- 
  full_data %>% 
  filter(.,
         !is.na(similarity_index_post2)) %>% 
  lm_lin(standardize(similarity_index_post2) ~ 
           treatment,
         covariates = ~  similarity_index_pre+
           block_id +male +
           na_smlr_blind_pre + na_smlr_imgrnt_pre +
           na_smlr_arb_pre + na_smlr_uo_pre,
         clusters = class_id,
         se_type = "stata",
         weights = assignment_prob_weights,
         data = .) 
o_sim_index_out2 <-
  tidy(o_sim_index2) %>% 
  filter(.,
         term == "treatment") %>% 
  mutate(.,
         Target = "Index",
         Outcome = "Group Similarity",
         Wave = "8-13 Weeks")



####Create Mechanism Plot-----

o_mech <- rbind(o_pt_index_out, o_pt_index_out2,
                o_sim_index_out, o_sim_index_out2,
                o_het_index_out, o_het_index_out2) %>% 
  mutate(.,
         se = paste0("(", round(std.error,2),")"),
         pese = paste(round(estimate,2), se))

o_mech$Outcome <- factor(o_mech$Outcome,
                         c("Group Heterogeneity", "Group Similarity", "Perspective Taking")) 

annotation_mech <- data.frame(
  x = c(-0.1,-0.1),
  y = c(3.15,2.85),
  label = c("1-2 Weeks", "8-13 Weeks")
)

ggplot(o_mech, aes(x=estimate, y=Outcome,
                   colour = forcats::fct_rev(Wave))) + 
  geom_point(position = position_dodge(width = 0.6))+
  geom_errorbar(aes(xmin=conf.low, xmax=conf.high), 
                width=.1,
                position = position_dodge(width = 0.6)) +
  geom_vline(xintercept = 0,
             color = "black",
             linetype = "dashed",
             size = 0.4) +
  xlim(-.18,.55)+
  geom_text(aes(label= pese, color = forcats::fct_rev(Wave)), 
            position = position_dodge(width = .6),
            hjust = -1.2, 
            size = 3,
            family = "Times",
            show.legend = FALSE) +
  geom_label(data=annotation_mech, aes(x=x, y=y, label=label),
             color=c("dodgerblue4", "firebrick3"), 
             size=2 , angle=45, fontface="bold")+
  labs(x = "",
       y = "",
       color = "")+
  scale_color_manual(values = c("firebrick3","dodgerblue4"))+
  # theme_tufte()+
  guides(color = guide_legend(reverse = TRUE))+
  theme(panel.spacing = unit(1.5, "lines"),
        text = element_text(size = 14, family = "Times"),
        # axis.text.y = element_text(size=14),
        legend.position = "none",
        panel.grid.major = element_blank(),
        panel.background = element_blank(),
        legend.key=element_blank())

##Table A7: Randomization Inference-----
###Thermometer First wave------
# Omit Outcome NAs and rename treatment and outcome
therm_p1_ri_data <-
  full_data %>% 
  filter(.,
         !is.na(therm_index_post)) %>% 
  mutate(.,
         Y = standardize(therm_index_post),
         Z = treatment)


# Declare assignment
declaration_therm_p1 <- 
  with(therm_p1_ri_data,{
    declare_ra(
      blocks = block_id,
      clusters = class_id)
  })

# Run RI
ri_therm_p1 <- 
  conduct_ri(
    Y ~ Z + therm_index_pre+
      block_id +male + 
      na_therm_arab_pre + na_therm_blind_pre + 
      na_therm_imgrnt_pre + na_therm_uo_pre,
    sharp_hypothesis = 0,
    IPW = T,
    IPW_weights = assignment_prob_weights,
    declaration = declaration_therm_p1,
    data = therm_p1_ri_data)
ri_therm_p1_out <- tidy(ri_therm_p1) %>% 
  mutate(.,
         term = "Thermometers W1")


###Thermometer Second wave------
# Omit Outcome NAs and rename treatment and outcome
therm_p2_ri_data <-
  full_data %>% 
  filter(.,
         !is.na(therm_index_post2)) %>% 
  mutate(.,
         Y = standardize(therm_index_post2),
         Z = treatment)


# Declare assignment
declaration_therm_p2 <- 
  with(therm_p2_ri_data,{
    declare_ra(
      blocks = block_id,
      clusters = class_id)
  })

# Run RI
ri_therm_p2 <- 
  conduct_ri(
    Y ~ Z + therm_index_pre+
      block_id +male + 
      na_therm_arab_pre + na_therm_blind_pre + 
      na_therm_imgrnt_pre + na_therm_uo_pre,
    sharp_hypothesis = 0,
    IPW = T,
    IPW_weights = assignment_prob_weights,
    declaration = declaration_therm_p2,
    data = therm_p2_ri_data)
ri_therm_p2_out <- tidy(ri_therm_p2) %>% 
  mutate(.,
         term = "Thermometers W2")



###Contact Intention First wave------
# Omit Outcome NAs and rename treatment and outcome
cntct_p1_ri_data <-
  full_data %>% 
  filter(.,
         !is.na(cntct_indx_post)) %>% 
  mutate(.,
         Y = standardize(cntct_indx_post),
         Z = treatment)


# Declare assignment
declaration_cntct_p1 <- 
  with(cntct_p1_ri_data,{
    declare_ra(
      blocks = block_id,
      clusters = class_id)
  })

# Run RI
ri_cntct_p1 <- 
  conduct_ri(
    Y ~ Z + cntct_indx_pre+
      block_id +male +
      na_cntct_1_blind_pre + na_cntct_2_blind_pre +
      na_cntct_3_blind_pre + na_cntct_4_blind_pre +
      na_cntct_1_arab_pre + na_cntct_2_arab_pre + 
      na_cntct_3_arab_pre + na_cntct_4_arab_pre +
      na_cntct_1_imgrnt_pre + na_cntct_2_imgrnt_pre +
      na_cntct_3_imgrnt_pre + na_cntct_4_imgrnt_pre +
      na_cntct_1_uo_pre + na_cntct_2_uo_pre +
      na_cntct_3_uo_pre + na_cntct_4_uo_pre,
    sharp_hypothesis = 0,
    IPW = T,
    IPW_weights = assignment_prob_weights,
    declaration = declaration_cntct_p1,
    data = cntct_p1_ri_data)
ri_cntct_p1_out <- tidy(ri_cntct_p1) %>% 
  mutate(.,
         term = "Contact W1")



###Contact Intention Second wave------
# Omit Outcome NAs and rename treatment and outcome
cntct_p2_ri_data <-
  full_data %>% 
  filter(.,
         !is.na(cntct_indx_post2)) %>% 
  mutate(.,
         Y = standardize(cntct_indx_post2),
         Z = treatment)


# Declare assignment
declaration_cntct_p2 <- 
  with(cntct_p2_ri_data,{
    declare_ra(
      blocks = block_id,
      clusters = class_id)
  })

# Run RI
ri_cntct_p2 <- 
  conduct_ri(
    Y ~ Z + cntct_indx_pre+
      block_id +male +
      na_cntct_1_blind_pre + na_cntct_2_blind_pre +
      na_cntct_3_blind_pre + na_cntct_4_blind_pre +
      na_cntct_1_arab_pre + na_cntct_2_arab_pre + 
      na_cntct_3_arab_pre + na_cntct_4_arab_pre +
      na_cntct_1_imgrnt_pre + na_cntct_2_imgrnt_pre +
      na_cntct_3_imgrnt_pre + na_cntct_4_imgrnt_pre +
      na_cntct_1_uo_pre + na_cntct_2_uo_pre +
      na_cntct_3_uo_pre + na_cntct_4_uo_pre,
    sharp_hypothesis = 0,
    IPW = T,
    IPW_weights = assignment_prob_weights,
    declaration = declaration_cntct_p2,
    data = cntct_p2_ri_data)
ri_cntct_p2_out <- tidy(ri_cntct_p2) %>% 
  mutate(.,
         term = "Contact W2")



###Diversity First wave------
# Omit Outcome NAs and rename treatment and outcome
div_p1_ri_data <-
  full_data %>% 
  filter(.,
         !is.na(div_scale_post)) %>% 
  mutate(.,
         Y = standardize(div_scale_post),
         Z = treatment)


# Declare assignment
declaration_div_p1 <- 
  with(div_p1_ri_data,{
    declare_ra(
      blocks = block_id,
      clusters = class_id)
  })

# Run RI
ri_div_p1 <- 
  conduct_ri(
    Y ~ Z + div_scale_pre+
      block_id +male + 
      na_div1_pre + na_div2_pre + na_div3_pre +
      na_div4_pre + na_div5_pre,
    sharp_hypothesis = 0,
    IPW = T,
    IPW_weights = assignment_prob_weights,
    declaration = declaration_div_p1,
    data = div_p1_ri_data)
ri_div_p1_out <- tidy(ri_div_p1) %>% 
  mutate(.,
         term = "Diversity W1")



###Diversity Second wave------
# Omit Outcome NAs and rename treatment and outcome
div_p2_ri_data <-
  full_data %>% 
  filter(.,
         !is.na(div_scale_post2)) %>% 
  mutate(.,
         Y = standardize(div_scale_post2),
         Z = treatment)


# Declare assignment
declaration_div_p2 <- 
  with(div_p2_ri_data,{
    declare_ra(
      blocks = block_id,
      clusters = class_id)
  })

# Run RI
ri_div_p2 <- 
  conduct_ri(
    Y ~ Z + div_scale_pre+
      block_id +male + 
      na_div1_pre + na_div2_pre + na_div3_pre +
      na_div4_pre + na_div5_pre,
    sharp_hypothesis = 0,
    IPW = T,
    IPW_weights = assignment_prob_weights,
    declaration = declaration_div_p2,
    data = div_p2_ri_data)
ri_div_p2_out <- tidy(ri_div_p2) %>% 
  mutate(.,
         term = "Diversity W2")

###Register First wave------
# Omit Outcome NAs and rename treatment and outcome
regcntct_p1_ri_data <-
  full_data %>% 
  filter(.,
         !is.na(contact_project)) %>% 
  mutate(.,
         Y = standardize(contact_project),
         Z = treatment)


# Declare assignment
declaration_regcntct_p1 <- 
  with(regcntct_p1_ri_data,{
    declare_ra(
      blocks = block_id,
      clusters = class_id)
  })

# Run RI
ri_regcntct_p1 <- 
  conduct_ri(
    Y ~ Z +  cntct_indx_pre+ therm_index_pre + div_scale_pre+
      block_id +male +
      na_cntct_1_blind_pre + na_cntct_2_blind_pre +
      na_cntct_3_blind_pre + na_cntct_4_blind_pre +
      na_cntct_1_arab_pre + na_cntct_2_arab_pre + 
      na_cntct_3_arab_pre + na_cntct_4_arab_pre +
      na_cntct_1_imgrnt_pre + na_cntct_2_imgrnt_pre +
      na_cntct_3_imgrnt_pre + na_cntct_4_imgrnt_pre +
      na_cntct_1_uo_pre + na_cntct_2_uo_pre +
      na_cntct_3_uo_pre + na_cntct_4_uo_pre+
      na_therm_arab_pre + na_therm_blind_pre + 
      na_therm_imgrnt_pre + na_therm_uo_pre+
      na_div1_pre + na_div2_pre + na_div3_pre +
      na_div4_pre + na_div5_pre,
    sharp_hypothesis = 0,
    IPW = T,
    IPW_weights = assignment_prob_weights,
    declaration = declaration_regcntct_p1,
    data = regcntct_p1_ri_data)
ri_regcntct_p1_out <- tidy(ri_regcntct_p1) %>% 
  mutate(.,
         term = "Register Contact W1")



###Register Second wave------
# Omit Outcome NAs and rename treatment and outcome
regcntct_p2_ri_data <-
  full_data %>% 
  filter(.,
         !is.na(contact_project_post2)) %>% 
  mutate(.,
         Y = standardize(contact_project_post2),
         Z = treatment)


# Declare assignment
declaration_regcntct_p2 <- 
  with(regcntct_p2_ri_data,{
    declare_ra(
      blocks = block_id,
      clusters = class_id)
  })

# Run RI
ri_regcntct_p2 <- 
  conduct_ri(
    Y ~ Z +  cntct_indx_pre+ therm_index_pre + div_scale_pre+
      block_id +male +
      na_cntct_1_blind_pre + na_cntct_2_blind_pre +
      na_cntct_3_blind_pre + na_cntct_4_blind_pre +
      na_cntct_1_arab_pre + na_cntct_2_arab_pre + 
      na_cntct_3_arab_pre + na_cntct_4_arab_pre +
      na_cntct_1_imgrnt_pre + na_cntct_2_imgrnt_pre +
      na_cntct_3_imgrnt_pre + na_cntct_4_imgrnt_pre +
      na_cntct_1_uo_pre + na_cntct_2_uo_pre +
      na_cntct_3_uo_pre + na_cntct_4_uo_pre+
      na_therm_arab_pre + na_therm_blind_pre + 
      na_therm_imgrnt_pre + na_therm_uo_pre+
      na_div1_pre + na_div2_pre + na_div3_pre +
      na_div4_pre + na_div5_pre,
    sharp_hypothesis = 0,
    IPW = T,
    IPW_weights = assignment_prob_weights,
    declaration = declaration_regcntct_p2,
    data = regcntct_p2_ri_data)
ri_regcntct_p2_out <- tidy(ri_regcntct_p2) %>% 
  mutate(.,
         term = "Register Contact W2")



###Bracelet Second wave------
# Omit Outcome NAs and rename treatment and outcome

brclt_ri_data <-
  full_data %>% 
  filter(.,
         school != "Democrati") %>% 
  mutate(.,
         Y = standardize(wrist_band),
         Z = treatment)


# Declare assignment
declaration_brclt <- 
  with(brclt_ri_data,{
    declare_ra(
      blocks = block_id,
      clusters = class_id)
  })

# Run RI
ri_brclt <- 
  conduct_ri(
    Y ~ Z + cntct_indx_pre+ therm_index_pre + div_scale_pre+ 
      block_id + male +
      na_cntct_1_blind_pre + na_cntct_2_blind_pre +
      na_cntct_3_blind_pre + na_cntct_4_blind_pre +
      na_cntct_1_arab_pre + na_cntct_2_arab_pre + 
      na_cntct_3_arab_pre + na_cntct_4_arab_pre +
      na_cntct_1_imgrnt_pre + na_cntct_2_imgrnt_pre +
      na_cntct_3_imgrnt_pre + na_cntct_4_imgrnt_pre +
      na_cntct_1_uo_pre + na_cntct_2_uo_pre +
      na_cntct_3_uo_pre + na_cntct_4_uo_pre+
      na_therm_arab_pre + na_therm_blind_pre + 
      na_therm_imgrnt_pre + na_therm_uo_pre+
      na_div1_pre + na_div2_pre + na_div3_pre +
      na_div4_pre + na_div5_pre,
    sharp_hypothesis = 0,
    IPW = T,
    IPW_weights = assignment_prob_weights,
    declaration = declaration_brclt,
    data = brclt_ri_data)
ri_brclt_out <- tidy(ri_brclt) %>% 
  mutate(.,
         term = "Bracelet W2")



###Perspective Taking First wave------
# Omit Outcome NAs and rename treatment and outcome
pt_p1_ri_data <-
  full_data %>% 
  filter(.,
         !is.na(pt_index_post)) %>% 
  mutate(.,
         Y = standardize(pt_index_post),
         Z = treatment)


# Declare assignment
declaration_pt_p1 <- 
  with(pt_p1_ri_data,{
    declare_ra(
      blocks = block_id,
      clusters = class_id)
  })

# Run RI
ri_pt_p1 <- 
  conduct_ri(
    Y ~ Z +  pt_index_pre+
      block_id +male +
      na_pt_arab_pre + na_pt_blind_pre +
      na_pt_imgrnt_pre + na_pt_uo_pre,
    sharp_hypothesis = 0,
    IPW = T,
    IPW_weights = assignment_prob_weights,
    declaration = declaration_pt_p1,
    data = pt_p1_ri_data)
ri_pt_p1_out <- tidy(ri_pt_p1) %>% 
  mutate(.,
         term = "Perspective Taking W1")



###Perspective Taking Second wave------
# Omit Outcome NAs and rename treatment and outcome
pt_p2_ri_data <-
  full_data %>% 
  filter(.,
         !is.na(pt_index_post2)) %>% 
  mutate(.,
         Y = standardize(pt_index_post2),
         Z = treatment)


# Declare assignment
declaration_pt_p2 <- 
  with(pt_p2_ri_data,{
    declare_ra(
      blocks = block_id,
      clusters = class_id)
  })

# Run RI
ri_pt_p2 <- 
  conduct_ri(
    Y ~ Z +  pt_index_pre+
      block_id +male +
      na_pt_arab_pre + na_pt_blind_pre +
      na_pt_imgrnt_pre + na_pt_uo_pre,
    sharp_hypothesis = 0,
    IPW = T,
    IPW_weights = assignment_prob_weights,
    declaration = declaration_pt_p2,
    data = pt_p2_ri_data)
ri_pt_p2_out <- tidy(ri_pt_p2) %>% 
  mutate(.,
         term = "Perspective Taking W2")



###Similarity First wave------
# Omit Outcome NAs and rename treatment and outcome
sim_p1_ri_data <-
  full_data %>% 
  filter(.,
         !is.na(similarity_index_post)) %>% 
  mutate(.,
         Y = standardize(similarity_index_post),
         Z = treatment)


# Declare assignment
declaration_sim_p1 <- 
  with(sim_p1_ri_data,{
    declare_ra(
      blocks = block_id,
      clusters = class_id)
  })

# Run RI
ri_sim_p1 <- 
  conduct_ri(
    Y ~ Z +  similarity_index_pre+
      block_id +male +
      na_smlr_blind_pre + na_smlr_imgrnt_pre +
      na_smlr_arb_pre + na_smlr_uo_pre,
    sharp_hypothesis = 0,
    IPW = T,
    IPW_weights = assignment_prob_weights,
    declaration = declaration_sim_p1,
    data = sim_p1_ri_data)
ri_sim_p1_out <- tidy(ri_sim_p1) %>% 
  mutate(.,
         term = "Similarity W1")



###Similarity second wave------
# Omit Outcome NAs and rename treatment and outcome
sim_p2_ri_data <-
  full_data %>% 
  filter(.,
         !is.na(similarity_index_post2)) %>% 
  mutate(.,
         Y = standardize(similarity_index_post2),
         Z = treatment)


# Declare assignment
declaration_sim_p1 <- 
  with(sim_p2_ri_data,{
    declare_ra(
      blocks = block_id,
      clusters = class_id)
  })

# Run RI
ri_sim_p2 <- 
  conduct_ri(
    Y ~ Z +  similarity_index_pre+
      block_id +male +
      na_smlr_blind_pre + na_smlr_imgrnt_pre +
      na_smlr_arb_pre + na_smlr_uo_pre,
    sharp_hypothesis = 0,
    IPW = T,
    IPW_weights = assignment_prob_weights,
    declaration = declaration_sim_p1,
    data = sim_p2_ri_data)
ri_sim_p2_out <- tidy(ri_sim_p2) %>% 
  mutate(.,
         term = "Similarity W2")




###Heterogeneity First wave------
# Omit Outcome NAs and rename treatment and outcome
het_p1_ri_data <-
  full_data %>% 
  filter(.,
         !is.na(het_index_post)) %>% 
  mutate(.,
         Y = standardize(het_index_post),
         Z = treatment)


# Declare assignment
declaration_het_p1 <- 
  with(het_p1_ri_data,{
    declare_ra(
      blocks = block_id,
      clusters = class_id)
  })

# Run RI
ri_het_p1 <- 
  conduct_ri(
    Y ~ Z + het_index_pre+
      block_id +male +
      na_htr_arab_pre + na_htr_blind_pre +
      na_htr_imgrnt_pre + na_htr_uo_pre,
    sharp_hypothesis = 0,
    IPW = T,
    IPW_weights = assignment_prob_weights,
    declaration = declaration_het_p1,
    data = het_p1_ri_data)
ri_het_p1_out <- tidy(ri_het_p1) %>% 
  mutate(.,
         term = "Heterogeneity W1")



###Heterogeneity Second wave------
# Omit Outcome NAs and rename treatment and outcome
het_p2_ri_data <-
  full_data %>% 
  filter(.,
         !is.na(het_index_post2)) %>% 
  mutate(.,
         Y = standardize(het_index_post2),
         Z = treatment)


# Declare assignment
declaration_het_p2 <- 
  with(het_p2_ri_data,{
    declare_ra(
      blocks = block_id,
      clusters = class_id)
  })

# Run RI
ri_het_p2 <- 
  conduct_ri(
    Y ~ Z + het_index_pre+
      block_id +male +
      na_htr_arab_pre + na_htr_blind_pre +
      na_htr_imgrnt_pre + na_htr_uo_pre,
    sharp_hypothesis = 0,
    IPW = T,
    IPW_weights = assignment_prob_weights,
    declaration = declaration_het_p2,
    data = het_p2_ri_data)
ri_het_p2_out <- tidy(ri_het_p2) %>% 
  mutate(.,
         term = "Heterogeneity W2")

####Create RI Table-------
ri_combined <-
  rbind(ri_therm_p1_out, ri_therm_p2_out, ri_cntct_p1_out, ri_cntct_p2_out,
        ri_div_p1_out, ri_div_p2_out, ri_regcntct_p1_out, ri_regcntct_p2_out,
        ri_brclt_out, ri_pt_p1_out, ri_pt_p2_out,ri_sim_p1_out, ri_sim_p2_out,
        ri_het_p1_out, ri_het_p2_out) %>% 
  rename(.,
         Term = term,
         Estimate = estimate) %>% 
  dplyr::select(.,
                Term:p.value) 
print(ri_combined)




## Figure A25: ATEs with wild cluster btstrp-------

####Thermometor------
therm_lm <- 
  full_data %>% 
  lm(standardize(therm_index_post) ~ 
       treatment + therm_index_pre+
       block_id +male + 
       na_therm_arab_pre + na_therm_blind_pre + 
       na_therm_imgrnt_pre + na_therm_uo_pre,
     weights = assignment_prob_weights,
     data = .) 

therm_bstrp <- 
  cluster.boot(therm_lm, full_data$class_id, 
               parallel = T, R = 1000, 
               boot_type = "wild")

therm_bstrp_out <- 
  coeftest(therm_lm, therm_bstrp) %>% 
  tidy(.,
       conf.int = T) %>% 
  mutate(.,
         Outcome = "Thermometers",
         Wave = "1-2 Weeks") %>% 
  filter(.,
         term == "treatment")


therm2_lm <- 
  full_data %>% 
  lm(standardize(therm_index_post2) ~ 
       treatment + therm_index_pre+
       block_id +male + 
       na_therm_arab_pre + na_therm_blind_pre + 
       na_therm_imgrnt_pre + na_therm_uo_pre,
     weights = assignment_prob_weights,
     data = .) 

therm2_bstrp <- 
  cluster.boot(therm2_lm, full_data$class_id, 
               parallel = T, R = 1000, 
               boot_type = "wild")

therm2_bstrp_out <- 
  coeftest(therm2_lm, therm2_bstrp) %>% 
  tidy(.,
       conf.int = T) %>% 
  mutate(.,
         Outcome = "Thermometers",
         Wave = "8-13 Weeks") %>% 
  filter(.,
         term == "treatment")


####Contact------
cntct_index_lm <- 
  full_data %>% 
  lm(standardize(cntct_indx_post) ~ 
       treatment + cntct_indx_pre+
       block_id +male +
       na_cntct_1_blind_pre + na_cntct_2_blind_pre +
       na_cntct_3_blind_pre + na_cntct_4_blind_pre +
       na_cntct_1_arab_pre + na_cntct_2_arab_pre + 
       na_cntct_3_arab_pre + na_cntct_4_arab_pre +
       na_cntct_1_imgrnt_pre + na_cntct_2_imgrnt_pre +
       na_cntct_3_imgrnt_pre + na_cntct_4_imgrnt_pre +
       na_cntct_1_uo_pre + na_cntct_2_uo_pre +
       na_cntct_3_uo_pre + na_cntct_4_uo_pre,
     weights = assignment_prob_weights,
     data = .) 


cntct_indx_bstrp <- 
  cluster.boot(cntct_index_lm, full_data$class_id, 
               parallel = T, R = 1000, 
               boot_type = "wild")

cntct_indx_bstrp_out <- 
  coeftest(cntct_index_lm, cntct_indx_bstrp) %>% 
  tidy(.,
       conf.int = T) %>% 
  mutate(.,
         Outcome = "Contact \nIntentions",
         Wave = "1-2 Weeks") %>% 
  filter(.,
         term == "treatment")

cntct2_index_lm <- 
  full_data %>% 
  lm(standardize(cntct_indx_post2) ~ 
       treatment + cntct_indx_pre+
       block_id +male +
       na_cntct_1_blind_pre + na_cntct_2_blind_pre +
       na_cntct_3_blind_pre + na_cntct_4_blind_pre +
       na_cntct_1_arab_pre + na_cntct_2_arab_pre + 
       na_cntct_3_arab_pre + na_cntct_4_arab_pre +
       na_cntct_1_imgrnt_pre + na_cntct_2_imgrnt_pre +
       na_cntct_3_imgrnt_pre + na_cntct_4_imgrnt_pre +
       na_cntct_1_uo_pre + na_cntct_2_uo_pre +
       na_cntct_3_uo_pre + na_cntct_4_uo_pre,
     weights = assignment_prob_weights,
     data = .) 


cntct2_indx_bstrp <- 
  cluster.boot(cntct2_index_lm, full_data$class_id, 
               parallel = T, R = 1000, 
               boot_type = "wild")

cntct2_indx_bstrp_out <- 
  coeftest(cntct2_index_lm, cntct2_indx_bstrp) %>% 
  tidy(.,
       conf.int = T) %>% 
  mutate(.,
         Outcome = "Contact \nIntentions",
         Wave = "8-13 Weeks") %>% 
  filter(.,
         term == "treatment")


####Contact Behavior------
cntct_bhvr_lm <- 
  full_data %>% 
  lm(standardize(contact_project) ~ 
       treatment +  cntct_indx_pre+ therm_index_pre + div_scale_pre+
       block_id +male +
       na_cntct_1_blind_pre + na_cntct_2_blind_pre +
       na_cntct_3_blind_pre + na_cntct_4_blind_pre +
       na_cntct_1_arab_pre + na_cntct_2_arab_pre + 
       na_cntct_3_arab_pre + na_cntct_4_arab_pre +
       na_cntct_1_imgrnt_pre + na_cntct_2_imgrnt_pre +
       na_cntct_3_imgrnt_pre + na_cntct_4_imgrnt_pre +
       na_cntct_1_uo_pre + na_cntct_2_uo_pre +
       na_cntct_3_uo_pre + na_cntct_4_uo_pre+
       na_therm_arab_pre + na_therm_blind_pre + 
       na_therm_imgrnt_pre + na_therm_uo_pre+
       na_div1_pre + na_div2_pre + na_div3_pre +
       na_div4_pre + na_div5_pre,
     weights = assignment_prob_weights,
     data = .) 

cntct_bhvr_bstrp <- 
  cluster.boot(cntct_bhvr_lm, full_data$class_id, 
               parallel = T, R = 1000, 
               boot_type = "wild")

cntct_bhvr_bstrp_out <- 
  coeftest(cntct_bhvr_lm, cntct_bhvr_bstrp) %>% 
  tidy(.,
       conf.int = T) %>% 
  mutate(.,
         Outcome = "Register \nContact",
         Wave = "1-2 Weeks") %>% 
  filter(.,
         term == "treatment")


cntct_bhvr2_lm <- 
  full_data %>% 
  lm(standardize(contact_project_post2) ~ 
       treatment +  cntct_indx_pre+ therm_index_pre + div_scale_pre+
       block_id +male +
       na_cntct_1_blind_pre + na_cntct_2_blind_pre +
       na_cntct_3_blind_pre + na_cntct_4_blind_pre +
       na_cntct_1_arab_pre + na_cntct_2_arab_pre + 
       na_cntct_3_arab_pre + na_cntct_4_arab_pre +
       na_cntct_1_imgrnt_pre + na_cntct_2_imgrnt_pre +
       na_cntct_3_imgrnt_pre + na_cntct_4_imgrnt_pre +
       na_cntct_1_uo_pre + na_cntct_2_uo_pre +
       na_cntct_3_uo_pre + na_cntct_4_uo_pre+
       na_therm_arab_pre + na_therm_blind_pre + 
       na_therm_imgrnt_pre + na_therm_uo_pre+
       na_div1_pre + na_div2_pre + na_div3_pre +
       na_div4_pre + na_div5_pre,
     weights = assignment_prob_weights,
     data = .) 

cntct_bhvr2_bstrp <- 
  cluster.boot(cntct_bhvr2_lm, full_data$class_id, 
               parallel = T, R = 1000, 
               boot_type = "wild")

cntct_bhvr2_bstrp_out <- 
  coeftest(cntct_bhvr2_lm, cntct_bhvr2_bstrp) %>% 
  tidy(.,
       conf.int = T) %>% 
  mutate(.,
         Outcome = "Register \nContact",
         Wave = "8-13 Weeks") %>% 
  filter(.,
         term == "treatment")

####Diversity------
div_index_lm <- 
  full_data %>% 
  lm(standardize(div_scale_post) ~ 
       treatment + div_scale_pre+
       block_id +male + 
       na_div1_pre + na_div2_pre + na_div3_pre +
       na_div4_pre + na_div5_pre,
     weights = assignment_prob_weights,
     data = .) 

div_bstrp <- 
  cluster.boot(div_index_lm, full_data$class_id, 
               parallel = T, R = 1000, 
               boot_type = "wild")

div_bstrp_out <- 
  coeftest(div_index_lm, div_bstrp) %>% 
  tidy(.,
       conf.int = T) %>% 
  mutate(.,
         Outcome = "Diversity \nAttitudes",
         Wave = "1-2 Weeks") %>% 
  filter(.,
         term == "treatment")


div_index2_lm <- 
  full_data %>% 
  lm(standardize(div_scale_post2) ~ 
       treatment + div_scale_pre+
       block_id +male + 
       na_div1_pre + na_div2_pre + na_div3_pre +
       na_div4_pre + na_div5_pre,
     weights = assignment_prob_weights,
     data = .) 

div2_bstrp <- 
  cluster.boot(div_index2_lm, full_data$class_id, 
               parallel = T, R = 1000, 
               boot_type = "wild")

div2_bstrp_out <- 
  coeftest(div_index2_lm, div2_bstrp) %>% 
  tidy(.,
       conf.int = T) %>% 
  mutate(.,
         Outcome = "Diversity \nAttitudes",
         Wave = "8-13 Weeks") %>% 
  filter(.,
         term == "treatment")


#### Diversity Behavior------

div_bhvr2_lm <- 
  full_data %>% 
  lm(standardize(wrist_band) ~ 
       treatment + cntct_indx_pre+ therm_index_pre + div_scale_pre+ 
       block_id + male +
       na_cntct_1_blind_pre + na_cntct_2_blind_pre +
       na_cntct_3_blind_pre + na_cntct_4_blind_pre +
       na_cntct_1_arab_pre + na_cntct_2_arab_pre + 
       na_cntct_3_arab_pre + na_cntct_4_arab_pre +
       na_cntct_1_imgrnt_pre + na_cntct_2_imgrnt_pre +
       na_cntct_3_imgrnt_pre + na_cntct_4_imgrnt_pre +
       na_cntct_1_uo_pre + na_cntct_2_uo_pre +
       na_cntct_3_uo_pre + na_cntct_4_uo_pre+
       na_therm_arab_pre + na_therm_blind_pre + 
       na_therm_imgrnt_pre + na_therm_uo_pre+
       na_div1_pre + na_div2_pre + na_div3_pre +
       na_div4_pre + na_div5_pre,
     weights = assignment_prob_weights,
     data = .) 


div_bhvr2_bstrp <- 
  cluster.boot(div_bhvr2_lm, full_data$class_id, 
               parallel = T, R = 1000, 
               boot_type = "wild")

div_bhvr2_bstrp_out <- 
  coeftest(div_bhvr2_lm, div_bhvr2_bstrp) %>% 
  tidy(.,
       conf.int = T) %>% 
  mutate(.,
         Outcome = "Diversity \nWristband",
         Wave = "8-13 Weeks") %>% 
  filter(.,
         term == "treatment")


##### Plot results with WCB------

summary_study2 <- rbind(
  therm_bstrp_out,therm2_bstrp_out,cntct_indx_bstrp_out,cntct2_indx_bstrp_out,
  cntct_bhvr_bstrp_out, cntct_bhvr2_bstrp_out,div_bstrp_out,div2_bstrp_out,
  div_bhvr2_bstrp_out) %>% 
  mutate(.,
         se = paste0("(", round(std.error,2),")"),
         pese = paste(round(estimate,2), se))

# set order
summary_study2$Outcome <- 
  factor(summary_study2$Outcome, 
         c("Diversity \nWristband", "Register \nContact", 
           "Diversity \nAttitudes",
           "Contact \nIntentions",
           "Thermometers"))

# Add annotation
annotation <- data.frame(
  x = c(-0.1,-0.1),
  y = c(5.2,4.8),
  label = c("1-2 Weeks", "8-13 Weeks")
)

ggplot(summary_study2, aes(x=estimate, y=Outcome,
                           colour = forcats::fct_rev(Wave))) + 
  geom_point(position = position_dodge(width = 0.6))+
  geom_errorbar(aes(xmin=conf.low, xmax=conf.high),
                position = position_dodge(width = 0.6),
                width=.1) +
  geom_vline(xintercept = 0,
             color = "black",
             linetype = "dashed",
             size = 0.4) +
  xlim(-0.18,.7)+
  #        geom_label()+
  geom_text(aes(label= pese, color = forcats::fct_rev(Wave)), 
            position = position_dodge(width = .6),
            family = "Times",
            hjust = -1.5, 
            size = 3,
            show.legend = FALSE) +
  geom_label(data=annotation, aes(x=x, y=y, label=label),
             color=c("dodgerblue4", "firebrick3"), 
             size=2 , angle=45, fontface="bold")+
  labs(x = "",
       y = "",
       color = "")+
  scale_color_manual(values = c("firebrick3", "dodgerblue4"))+
  # theme_tufte()+
  guides(color = guide_legend(reverse = TRUE))+
  theme(panel.spacing = unit(1.5, "lines"),
        text = element_text(size = 15, family = "Times"),
        # axis.text.y = element_text(size=14),
        legend.position = "none",
        panel.grid.major = element_blank(),
        panel.background = element_blank(),
        legend.key=element_blank())


## Figure A26: Mechanisms with wild cluster btstrp-----

###Perspective Taking-----
pt_index_lm <- 
  full_data %>% 
  lm(standardize(pt_index_post) ~ 
       treatment + pt_index_pre+
       block_id +male +
       na_pt_arab_pre + na_pt_blind_pre +
       na_pt_imgrnt_pre + na_pt_uo_pre,
     weights = assignment_prob_weights,
     data = .) 


pt_bstrp <- 
  cluster.boot(pt_index_lm, full_data$class_id, 
               parallel = T, R = 1000, 
               boot_type = "wild")

pt_bstrp_out <- 
  coeftest(pt_index_lm, pt_bstrp) %>% 
  tidy(.,
       conf.int = T) %>% 
  mutate(.,
         Outcome = "Perspective Taking",
         Wave = "1-2 Weeks") %>% 
  filter(.,
         term == "treatment")



pt_index2_lm <- 
  full_data %>% 
  lm(standardize(pt_index_post2) ~ 
       treatment + pt_index_pre+
       block_id +male +
       na_pt_arab_pre + na_pt_blind_pre +
       na_pt_imgrnt_pre + na_pt_uo_pre,
     weights = assignment_prob_weights,
     data = .) 

pt2_bstrp <- 
  cluster.boot(pt_index2_lm, full_data$class_id, 
               parallel = T, R = 1000, 
               boot_type = "wild")

pt2_bstrp_out <- 
  coeftest(pt_index2_lm, pt2_bstrp) %>% 
  tidy(.,
       conf.int = T) %>% 
  mutate(.,
         Outcome = "Perspective Taking",
         Wave = "8-13 Weeks") %>% 
  filter(.,
         term == "treatment")


### Heterogeneity------

het_index_lm <- 
  full_data %>% 
  lm(standardize(het_index_post) ~ 
       treatment + het_index_pre+
       block_id +male +
       na_htr_arab_pre + na_htr_blind_pre +
       na_htr_imgrnt_pre + na_htr_uo_pre,
     weights = assignment_prob_weights,
     data = .) 

het_bstrp <- 
  cluster.boot(het_index_lm, full_data$class_id, 
               parallel = T, R = 1000, 
               boot_type = "wild")

het_bstrp_out <- 
  coeftest(het_index_lm, het_bstrp) %>% 
  tidy(.,
       conf.int = T) %>% 
  mutate(.,
         Outcome = "Group Heterogeneity",
         Wave = "1-2 Weeks") %>% 
  filter(.,
         term == "treatment")



het_index2_lm <- 
  full_data %>% 
  lm(standardize(het_index_post2) ~ 
       treatment + het_index_pre+
       block_id +male +
       na_htr_arab_pre + na_htr_blind_pre +
       na_htr_imgrnt_pre + na_htr_uo_pre,
     weights = assignment_prob_weights,
     data = .) 

het2_bstrp <- 
  cluster.boot(het_index2_lm, full_data$class_id, 
               parallel = T, R = 1000, 
               boot_type = "wild")

het2_bstrp_out <- 
  coeftest(het_index2_lm, het2_bstrp) %>% 
  tidy(.,
       conf.int = T) %>% 
  mutate(.,
         Outcome = "Group Heterogeneity",
         Wave = "8-13 Weeks") %>% 
  filter(.,
         term == "treatment")



### Similarity------

sim_index_lm <- 
  full_data %>% 
  lm(standardize(similarity_index_post) ~ 
       treatment +similarity_index_pre+
       block_id +male +
       na_smlr_blind_pre + na_smlr_imgrnt_pre +
       na_smlr_arb_pre + na_smlr_uo_pre,
     weights = assignment_prob_weights,
     data = .) 

sim_bstrp <- 
  cluster.boot(sim_index_lm, full_data$class_id, 
               parallel = T, R = 1000, 
               boot_type = "wild")

sim_bstrp_out <- 
  coeftest(sim_index_lm, sim_bstrp) %>% 
  tidy(.,
       conf.int = T) %>% 
  mutate(.,
         Outcome = "Group Similarity",
         Wave = "1-2 Weeks") %>% 
  filter(.,
         term == "treatment")



sim_index2_lm <- 
  full_data %>% 
  lm(standardize(similarity_index_post2) ~ 
       treatment +similarity_index_pre+
       block_id +male +
       na_smlr_blind_pre + na_smlr_imgrnt_pre +
       na_smlr_arb_pre + na_smlr_uo_pre,
     weights = assignment_prob_weights,
     data = .) 

sim2_bstrp <- 
  cluster.boot(sim_index2_lm, full_data$class_id, 
               parallel = T, R = 1000, 
               boot_type = "wild")

sim2_bstrp_out <- 
  coeftest(sim_index2_lm, sim2_bstrp) %>% 
  tidy(.,
       conf.int = T) %>% 
  mutate(.,
         Outcome = "Group Similarity",
         Wave = "8-13 Weeks") %>% 
  filter(.,
         term == "treatment")



###Create Mechanism Plot-----

mech <- rbind(pt_bstrp_out,pt2_bstrp_out,
              het_bstrp_out,het2_bstrp_out,
              sim_bstrp_out,sim2_bstrp_out) %>% 
  mutate(.,
         se = paste0("(", round(std.error,2),")"),
         pese = paste(round(estimate,2), se))

mech$Outcome <- factor(mech$Outcome,
                       c("Group Heterogeneity", "Group Similarity", "Perspective Taking")) 

annotation_mech <- data.frame(
  x = c(-0.1,-0.1),
  y = c(3.15,2.85),
  label = c("1-2 Weeks", "8-13 Weeks")
)

ggplot(mech, aes(x=estimate, y=Outcome,
                 colour = forcats::fct_rev(Wave))) + 
  geom_point(position = position_dodge(width = 0.6))+
  geom_errorbar(aes(xmin=conf.low, xmax=conf.high), 
                width=.1,
                position = position_dodge(width = 0.6)) +
  geom_vline(xintercept = 0,
             color = "black",
             linetype = "dashed",
             size = 0.4) +
  xlim(-.18,.55)+
  geom_text(aes(label= pese, color = forcats::fct_rev(Wave)), 
            position = position_dodge(width = .6),
            hjust = -1.5, 
            size = 3,
            family = "Times",
            show.legend = FALSE) +
  geom_label(data=annotation_mech, aes(x=x, y=y, label=label),
             color=c("dodgerblue4", "firebrick3"), 
             size=2 , angle=45, fontface="bold")+
  labs(x = "",
       y = "",
       color = "")+
  scale_color_manual(values = c("firebrick3","dodgerblue4"))+
  # theme_tufte()+
  guides(color = guide_legend(reverse = TRUE))+
  theme(panel.spacing = unit(1.5, "lines"),
        text = element_text(size = 14, family = "Times"),
        # axis.text.y = element_text(size=14),
        legend.position = "none",
        panel.grid.major = element_blank(),
        panel.background = element_blank(),
        legend.key=element_blank())


## Figure A27 panels a/b/c/d: External Validity Bias----
###Thermometer-----
full_data_therm_exr <-
  full_data %>% 
  filter(.,
         !is.na(therm_index_post))

covariates <- c("male", "block_id", "age",
                "therm_index_pre", "div_scale_pre", "cntct_indx_pre",
                "pt_index_pre", "het_index_pre", "similarity_index_pre") 
for_sate_therm <- as.formula(paste0("therm_index_post ~ treatment + ", 
                                    paste(covariates, collapse = "+")))

lm_sate_therm <- lm_robust(for_sate_therm, 
                           clusters = class_id,
                           se_type = "stata",
                           data = full_data_therm_exr)
sate_est_therm <- summary(lm_sate_therm)$coef["treatment", c(1,2)]
therm_exr <- exr(outcome = "therm_index_post", 
                 treatment = "treatment", 
                 covariates = c("male", "block_id", "age",
                                "therm_index_pre", "div_scale_pre", "cntct_indx_pre",
                                "pt_index_pre", "het_index_pre", "similarity_index_pre"),
                 data = full_data_therm_exr,
                 uncertainty = TRUE,
                 numCores = 4,
                 clusters = full_data_therm_exr$class_id,
                 boot = T,
                 sate_estimate = sate_est_therm)
plot(therm_exr)
dev.off() 

#### Contact-----
full_data_cnt_exr <-
  full_data %>% 
  filter(.,
         !is.na(cntct_indx_post))

for_sate_cnt <- as.formula(paste0("cntct_indx_post ~ treatment + ", 
                                  paste(covariates, collapse = "+")))

lm_sate_cnt <- lm_robust(for_sate_cnt, 
                         clusters = class_id,
                         se_type = "stata",
                         data = full_data_cnt_exr)
sate_est_cnt <- summary(lm_sate_cnt)$coef["treatment", c(1,2)]
cnt_exr <- exr(outcome = "cntct_indx_post", 
               treatment = "treatment", 
               covariates = c("male", "block_id", "age",
                              "therm_index_pre", "div_scale_pre", "cntct_indx_pre",
                              "pt_index_pre", "het_index_pre", "similarity_index_pre"),
               data = full_data_cnt_exr,
               uncertainty = TRUE,
               numCores = 4,
               clusters = full_data_cnt_exr$class_id,
               boot = T,
               sate_estimate = sate_est_cnt)
plot(cnt_exr)
dev.off() 



#### Diveristy-----
full_data_div_exr <-
  full_data %>% 
  filter(.,
         !is.na(div_scale_post))

for_sate_div <- as.formula(paste0("div_scale_post ~ treatment + ", 
                                  paste(covariates, collapse = "+")))

lm_sate_div <- lm_robust(for_sate_div, 
                         clusters = class_id,
                         se_type = "stata",
                         data = full_data_div_exr)
sate_est_div <- summary(lm_sate_div)$coef["treatment", c(1,2)]
div_exr <- exr(outcome = "div_scale_post", 
               treatment = "treatment", 
               covariates = c("male", "block_id", "age",
                              "therm_index_pre", "div_scale_pre", "cntct_indx_pre",
                              "pt_index_pre", "het_index_pre", "similarity_index_pre"),
               data = full_data_div_exr,
               uncertainty = TRUE,
               numCores = 4,
               clusters = full_data_div_exr$class_id,
               boot = T,
               sate_estimate = sate_est_div)
plot(div_exr)
dev.off() 



#### Diveristy Wristband-----
full_data_wb_exr <-
  full_data %>% 
  filter(.,
         !is.na(wrist_band))

for_sate_wb <- as.formula(paste0("wrist_band ~ treatment + ", 
                                 paste(covariates, collapse = "+")))

lm_sate_wb <- lm_robust(for_sate_wb, 
                        clusters = class_id,
                        se_type = "stata",
                        data = full_data_wb_exr)
sate_est_wb <- summary(lm_sate_wb)$coef["treatment", c(1,2)]
wb_exr <- exr(outcome = "wrist_band", 
              treatment = "treatment", 
              covariates = c("male", "block_id", "age",
                             "therm_index_pre", "div_scale_pre", "cntct_indx_pre",
                             "pt_index_pre", "het_index_pre", "similarity_index_pre"),
              data = full_data_wb_exr,
              uncertainty = TRUE,
              numCores = 4,
              clusters = full_data_wb_exr$class_id,
              boot = T,
              sate_estimate = sate_est_wb)
plot(wb_exr)
dev.off() 

